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Abstract 


A  high-order  stabilization  filter  was  formerly  designed  to  stabilize 
an-unstable  pitch  control  system  of  a. terminal  homing  missile  system.  In 
this  report,  a  new  dominant-data  matching  method  is  presented  to  redesign 
the  high  order  stabilization  filter.  Using  this  new  method  several  re¬ 
duced  order  filters  are  obtained.  As  a  result,  the  cost  of  implementation 
is  reduced  and  the  reliability  is  increased.  An  algebraic  method  is  also 
applied  to  redesign  the  stabilization  filter  so  that  the  performance  of 
the  redesigned  pitch  control  system  is  improved.  In  addition,  the  pro¬ 
posed  dominant-data  matching  method  can  be  applied  to  determine  a  reduced 
order  model  of  a  high  order  system.  Unlike  the  reduced  order  models  ob¬ 
tained  by  most  existing  model  reduction  methods,  the  reduced  order  model 
mentioned  above  has  the  exact  assigned  frequency- domain  specifications  of 
the  original  system.  The  dominant-data  matching  method  can  also  be  applied 
to  identify  any  practical  system.  , 


ABSTRACT 


A  high-order  stabilization  filter  was  formerly  designed  to 
stabilize  an  unstable  pitch  control  system  of  a  terminal  homing  missile 
system.  In  this  report,  a  new  dominant-data  matching  method  is  pre¬ 
sented  to  redesign  the  high  order  stabilization  filter.  Using  this  new 
method  several  reduced  order  filters  are  obtained.  As  a  result,  the 
cost  of  implementation  is  reduced  and  the  reliability  is  increased.  An 
algebraic  method  is  also  applied  to  redesign  the  stabilization  filter 
so  that  the  performance  of  the  redesigned  pitch  control  system  is  im¬ 
proved.  In  addition,  the  proposed  dominant-data  matching  method  can  be 
applied  to  determine  a  reduced  order  model  of  a  high  order  system. 
Unlike  the  reduced  order  models  obtained  by  most  existing  model  reduc- 
tion  methods,  the  reduced  order  model  mentioned  above  has  the  exact 
assigned  frequency-domain  specifications  of  the  original  system.  The 
dominant-data  matching  method  can  also  be  applied  to  identify  any 
practical  system. 
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CHAPTER  I 


INTRODUCTION 

This  report  deals  with  the  simplification  and  realization  of  a 
stabilization  filter  designed  to  stabilize  the  pitch  control  system  of 
an  unstable  semi-active  terminal  homing  missile  system  [1].  The  block 
diagram  of  the  existing  stabilized  system  is  shown  in  Fig.  1. 


Figure  1.  The  Block  Diagram  of  the  Existing  Control  System 


The  overall  transfer  function  of  the  existing  system  shown  in  Fig.  1 
is  given  by 
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=  The  transfer  function  of  the  actuator  and  the  air  frame  dynamics 


of  the  missile  system. 
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=  The  open  loop  transfer  function  of  the  existing  stabilized  system. 


H  (s)  =  Transfer  function  of  the  gyro, 

g 


=  1,  as  the  rate  gyro  is  not  present  in  the  system. 


After  substituting  H^(s) 


1  and  Eqn.  (1.2)  and  (1.3)  into  Eqn.  (1.1) 


it  becomes 
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From  Fig.  1  as  well  as  from  Eqn.  (1.2)  it  can  be  noticed  that  the  exist¬ 
ing  stabilization  filter  jrstaj,(s)  is  a  fourth-order  series  compensator 

with  two  pairs  of  complex  poles.  F  .  (s)  is  not  a  positive  real  function 

stab 

and  hence  cannot  be  synthesized  with  passive  elements.  The  objective  of 


this  report  is  to  develop  computer-aided  design  methods  for  redesigning 
the  stabilization  filter  in  a  simpler  form  so  that  the  cost  of  imple¬ 
mentation  of  the  compensator  can  be  reduced  and  at  the  same  time  the 
performance  of  the  redesigned  pitch  control  system  remains  the  same  as 
that  of  the  existing  pitch  control  system. 

Nyquist  plots  of  G£(s)  and  Gq(s)  are  shown  in  Fig.  2.  The 
dominant  frequency-response  data  of  G£(s)  are  given  below: 


i)  The  real  and  imaginary  parts  of  Ge(s)  at  s  =  jw  =  jO  are 

Re  [Gg(jO)]  =  -2.103817  and  1m  [Gg(j0)]  =  ®  (2.1) 

or  Tg(j0)  =  1 

ii)  The  gain  margin  of  this  system  Ge(j<o)  is 
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where  the  phase-crossover  frequency  w  is  given  by 
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Re [C  (jw  )]  =  -1.507944  (2.4) 
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iii)  The  phase  margin  4>em  of  the  system  G^Cjw)  is 
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where  the  gain  cross-over  frequency  w  is  given  by  -3.2 
rad/sec  so  that 


=  1 


(2.7) 


The  equivalent  real  and  imaginary  parts  of  G  (jw)  at  u  *  uec 


3.2  rad/sec.  are 
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The  frequency  response  data  at  w  ■  0  in  (2.1)  indirectly  indi¬ 
cates  the  steady-state  value  of  the  unit  step  response  of  the  system 

T  (s)  .  The  data  at  uj  =  u>  and  to  =  u  in  Eqn.  (2)  represent  two  con- 
e  eir  ec 

trol  specifications  [2]:  gain  margin  and  phase  margin.  These  control 
specifications  characterize  the  relative  stability  and  the  transient  re¬ 
sponse  of  the  existing  stabilized  system.  The  dominant  frequency  response 

data  of  G.(s),  F  ,  (s)  etc.  are  listed  below: 
u  stab 

i)  The  real  and  imaginary  parts  of  Gp(jw)  at  u  =  0  are 
Re [GQ (jO) ] =  -1.304841  and  Im(tyjO) J  =  » 
ii)  The  phase  margin  4>Qm  of  the  original  system  G0(jio)  is 
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(3.1) 


*  ••  V 


•  7 


K0m 


=  180°  +  /G0(jw0c)  =  -5.58* 


(3.2) 


where  the  gain  crossover  frequency  is  given  by 


-  1.6  rad/sec.  so  that  |G.(ju>_  )  I  =  1.  (3.3) 

Uc  0  0c 


Other  frequency  response  data  at  (±>e^  =  1.9  rad/sec.  and  =  3.2 


rad /sec.  are 
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iv)  Re[Gn(ju)  )]  =  -0.6181657 
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The  dominant  frequency  response  data  of  the  stabilization  filter 

F  ,(s)  are 
stab 
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or 

lFstab«“J  ‘  1-614464333  and  /Fstab(3“ec>  ~  7.293349493’ 

at  w  =  3.2  rad/sec.  (4.5) 

ec 

Now,  analyzing  the  data  we  have  from  Eqn.  (1.3)  and  (1.4),  it  is  clear 

that  G,,(s)  and  G  (s)  are  non-minimum  phase  functions  and  thev  are  un- 
u  e 

stable  because  of  the  pole  s  =  2.821  which  is  in  the  right  half  plane 

of  the  s-plane.  Referring  to  the  Nyquist  plots  in  Fig.  2,  and  according  to 

Nyquist  stability  criterion  the  original  missile  system  (without  F  ,(s) ) 

stab 

is  unstable  whereas  the  existing  stabilized  system  is  asymptotically 
stable.  However,  due  to  the  small  positive  phase  margin  given  in  Eqn. 
(2.5),  the  time  response  of  the  existing  stabilized  system  is  quite  os¬ 
cillatory. 

To  redesign  the  pitch  control  system  or  the  stabilization  filter 
s>  that  the  cost  of  implementation  is  reduced  and  the  flight  control 
performance  of  the  missile  system  is  improved,  two  computer-aided  methods 
have  been  developed.  These  two  proposed  methods  will  be  discussed  in 
this  report  Step  by  step.  In  Chapter  II  a  transfer  function  (called  a 
standard  transfer  function  Tr(s))  is  obtained  by  using  a  dominant-data 
matching  method.  Tr(s)  matches  the  assigned  specifications  given  in 
Eqn.  (2).  Therefore,  the  standard  transfer  function  Tr(s)  mentioned 
above  is  a  reduced -order  model  of  the  existing  stabilized  system  T^(s) 
in  Eqn.  (1.5).  The  unit  step  response  curves  of  Tg(s)  and  Tf(s) 

will  be  compared.  This  comparison  will  also  verify  that  the  data  in  Eqn. 
(2)  are  dominant  ones. 

To  solve  the  nonlinear  equations  obtained  in  Chapter  II  four 
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different  methods  of  finding  initial  guesses  are  discussed  in  Chapter 
III. 

In  Chapter  IV  two  reduced  order  models  of  the  stabilization  fil¬ 
ter  Fstaj,(s)  are  obtained.  One  of  these  two  is  obtained  by  using  the 
dominant-data  matching  method  and  the  other  by  using  a  similar  approach 
to  fit  a  low  order  model  that  satisfies  the  specifications  shown  in  Eqn. 

(4). 

Chapter  V  consists  of  two  parts,  in  the  first  part  the  dominant- 
data  matching  method  is  used  to  obtain  an  unstable  reduced  order  model 
of  the  original  high-order  unstable  system  G^(s)  shown  in  Eqn.  (1.3). 

This  is  done  just  to  simplify  the  design  procedure.  In  the  second  part 
of  Chapter  V  the  algebraic  method  of  Shieh  [3]  and  Chen  [4]  is  applied  to  re¬ 
design  the  pitch  control  system.  This  is  done  by  designing  a  series 
filter  in  the  feed  forward  path  and  a  parallel  filter  in  the  feedback 
path.  Thus,  the  advantages  of  feedback  control  structure  have  been  fully 


utilized. 


CHAPTER  II 


THE  DOMINANT-DATA  MATCHING  METHOD 

The  design  goals  of  a  control  system  are  often  characterized  by 
a  set  of  control  specifications.  These  specifications  can  be  classi¬ 
fied  as  i)  time-domain  specifications  such  as  rise  time,  setting  time, 
ii)  frequency  domain  specifications  such  as  phase  margin,  gain  margin 
and  iii)  complex  domain  specifications  such  as  damping  ratio,  and  natural 
angular  frequency.  Empirical  rules  that  link  the  specifications  in  the 
time,  frequency,  and  complex  domains  have  been  developed  by  Truxal 
[5],  Toro  and  Parker  [6],  Axelby  [7]  and  Seshadri  et.  al.  [8].  From 
these  results,  it  is  observed  that  most  time-domain  specifications  and 
complex-domain  specifications  can  be  approximately  converted  to  fre¬ 
quency-domain  specifications.  Some  of  these  frequency-domain  specifi¬ 
cations  are  phase  margin  (4>m)  ,  gain  margin  (G^) ,  maximum  value  of  the 
closed-loop  frequency  response  (M^),  phase  crossover  frequency  (w^), 

gain-crossover  frequency  (u  ),  peak  value  frequency  (w  ),  the  bandwidth 

c  p 

(w^) ,  and  the  velocity  error  constant  Other  important  frequency  re¬ 

sponse  data  are: 

(1)  The  real  part  and  imaginary  part  of  the  closed-loop  function 
T(s)  as  well  as  the  corresponding  open-loop  function  G(s)  at 
s  =  joj  =  jO, 

(2)  the  real  part  of  the  open-loop  transfer  function  G(ju)  at  the 
phase  crossover  frequency  which  has  been  used  to  define  the 
gain  margin  (G  ) , 


(3)  the  corner  frequencies  in  the  Bode  plot  of  G(jw)  in  the  regions 
of  to  =  wc  where  20  logjGfjw^)!  »  +  15  db  and  u  =  u  ^  where 
20  log|C(jwc2)|  =  -15  db. 

Chen  [9]  has  shown  empirically  that  the  open-loop  poles  and  zeros  of  a 
system  can  be  approximated  by  retaining  the  Bode  plot  in  the  regions  of 
the  ±  15  db  boundaries. 

The  data  at  co  =  0  often  indicate  the  final  value  and  the  type 
of  the  system.  More  specifically,  the  value  of  T(j0)  or  real  part  of  G(j0)  indi 
cates  the  final  value  of  the  system,  while  the  imaginary  part  of  G(j0)  indi¬ 
cates  the  type  of  the  system.  For  example,  for  a  type  'O'  system,  the 
imaginary  part  of  G(j0)  is  0,  and  for  any  system  other  than  type  'O', 
for  example,  say  type  '1',  it  is  infinity. 

Depending  on  the  problem  one  has,  one  can  use  any  one  or  a  com¬ 
bination  of  the  time  domain,  frequency  domain  and  complex-domain  speci¬ 
fications.  However,  in  this  case  the  original  pitch  control  system  that 
is  available  is  a  high  order  transfer  function  with  large  coefficients 
Eqn.(1.5).  As  a  result  the  time  response  curve  and  the  corresponding 
time  domain  specifications  of  this  practical  system  T£(s)  are  difficult 
to  obtain.  On  the  other  hand,  with  the  help  of  a  digital  computer  the 
frequency  response  curve  and  hence  the  frequency  domain  specifications 
of  Te(s)  can  be  easily  determined.  Therefore,  a  frequency  domain  approach 
or  a  dominant  data  matching  method  is  proposed  to  construct  a  transfer 
function  Tr(s),  a  reduced  order  model  of  Tg(s),  and  to  redesign  the 
pitch  control  system.  Several  methods  have  been  already  proposed  (10, 

11,12]  to  obtain  reduced  order  models  by  considering  frequency  domain 
specifications.  However,  the  only  reduced  order  models  that  satisfy  the 
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assigned  specifications  exactly  are  the  ones  obtained  by  the  proposed 


method. 


From  the  rules  of  thumb  it  is  observed  that  the  gain  margin, 


the  phase  margin,  the  gain  cross-over  frequency  and  the  phase  cross¬ 
over  frequency  are  the  most  important  frequency  response  data.  These 
data  are  called  the  dominant  frequency  response  data.  Another  impor¬ 
tant  frequency  response  data  is  the  steady  state  value  of  a  closed-loop 
system,  which  is  indirectly  given  by  the  value  of  the  real  part  of  the 
open  loop  transfer  function  G ( j to)  at  u  *  0.  These  dominant  frequency 
response  data  may  be  considered  as  the  design  goal.  Let  us  assume  that 

the  desired  reduced  order  model  of  T  (s)  which  may  be  called  the  standard 

e 

model  of  T  (s)  is 
e 


b+b  s+b  s 

y»>  =  ---—2 — 3 

a^+a^s+a^s  +a^s 


It  is  required  that  T^(s)  satisfies  all  the  conditions  given 
in  Eqn.  (2) . 

From  the  conditions  in  (2.1),  it  may  be  observed  that  the  sys¬ 
tem  T^(s)  is  a  type  1  system.  Therefore  b^  =  a^.  Also,  to  simplify  the 


equation  we  let  a^  =  1.  Thus,  we  have 


a.+b  s+b  s  G  (s) 

Tr(s)  =  •  r  273  -  w-go 

a^i-a^si^s  +s  r 


where  the  open-loop  transfer  function  G^(s)  is  given  by 

2  2 
a  +b^s+b2S  K[1+CjS+C2S  ] 

G  (s)  - - “  27 

s[(a  -bjH^-b^s+s  ]  s[l+d^s+d2s  ] 


(5.1) 


(5.2) 
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where 


K  = 


arbi  * 


a2~b2  1 

d!  -I^bT  •  d2=dr 


*1  “1 


1  1 


The  unknown  coefficients  and  are  to  be  determined  by  us¬ 

ing  the  conditions  in  Eqn.  (2).  Following  the  basic  definitions  and 
substituting  the  assigned  data  in  Eqn.  (2)  yields  a  set  of  nonlinear 
equations  f ^ (a^ ,a^ ,a^ ,b^ ,b2 )  =  0  for  i  =  1,2,. ..,5  as  follows: 

K[  1+juic  +(jus)2c  ] 

i)  G  (jo»)  =  - 2 - 

jw[  l+ja)d^+(  joi)  d2] 

=  [l+ju(c  -d1)+(jw)^(  )  +  ...  ] 

jw  1  1 

a  >  ll+j“(crdi)] 

=  KCc^)  -  j  l 

i>o  Gr(ju)  =  K(crdi) ' 


K*Co  Gr<^>!  “  K(cl'dl)  =7^1 


Eqn.  (2.1)  gives  RelG^CjO)]  =  -2.1 


arY 


(6.0) 


or 


or 


or 


(bi  . 

arbi  ao  arbi 


i 


.  a0(a2~b2) 
al"bl  (a1~b1)2 


-2.1 


-2.1 


Ya0,al’a2’bl*Y  ’  Yal”Y'Ya2“Y+2,1^al~Y  »  0 

(6.1) 


-1.5,  at  u> 


1.9  rad/sec 


The  data  in  (2.2), 


or  Re[G  (jw  )] 
r  eTT 


eu 


gives 


Re[Gr(jw)J 


Ref 


(aQ-u>  b2)+jujb1 


-] 


UJ—tiJ 


eir 


-u )2  (a2-b2)+jw(a^-b^-a)^) 

w=u  =1.9 
ew 


-1.5 


or 

(a n-w2  b»)(a  -b-)+w2  b  (a  -b.-u)2  ) 
0  er  2  2  2  eir  1  1  1 

4  ,  .  .2.  2  ,  .  2.2 

u>  (a„-b„)  +w  (a  -b  -a)  ) 

eir  2  2  eir  l  l  eTT 


1.5 


w  =1.9 

eTi 


or 


f2^aO,al,a2,bl,b2^  =  ^a2-b2^a0_3'61b2^-bl^ai“bl”^'^1^ 
-  1.5[3.61(a2-b2)2+(a1-b1-3.61)2]  = 


The  condition  in  (2.3),  or  =  -180°  where  w 

rad /sec,  gives 


eTi 


tan 


-1 


a)  b, 
eir  1 

VWeTTb2 


180°  +  tan 


■1  ^eir (  3  l~b  l~tl>eir  > 

“e^VV 


-180“ 


or 


tan 


1 . 9b  j  aj-bj-3.61 

-i  v3-61l,2 + 


1  -  -r 


1.9b.,  (a1~b1~3. 61) 
l79(a0-3.61b2)(a2-b2) 


=  0° 


or 


f3^a0,a] ’a2,bl,b2^  =  3.61b1(a2-b2)+(aQ-3.61b2)  (aj^-bj^-3 


The  data  in  (2.6)  or 


em 


180° 


5.7787°,  yields 


w  =3.2  rad/sec 
ec 


0  (6.2) 

=  1.9 


.61)  =  0 
(6.3) 
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180°  +  tan 


-1 


3.2b, 


a0-10.24b2 


180*  +  tan 


-i  arbr10,24 

3.2(a2-b2) 


5 . 7787* 


or 


3.2b, 


tan 


_j  a0~10.24b2 


a1-b1-10.24 

3.2(a2-b2) 


3.2b1(a1-b1-10.24) 


=  5.7787° 


1  - 


3.2(a2-b2)(a0-10.24b2) 


or 


10. 24b1(a2-b2)+(a()-10. 2482X3^^10. 24) 
y.2(a2-b2)(a0-10.24b2)-3.2b1(a1-b1-10.24) 


0.10120072 


or 


f4(a0,a1,a2,b1,b2)  =  10. 24bx (a2-b2)+(a0-10. 24b2)  (a^-10. 24) 

-  0.  3238423014  [  Uj-b^  (aQ-10. 24b2>-b1(a;L-b1 

-  10.24)]  =  0  (6.4) 

v)  The  condition  in  (2.7)  or 


|G  (iio  )|  =  1  where  u  =  3.2  rad/sec,  gives 
1  r  ec  1  ec 

a  -10. 24b„+j  3.2b, 

1 _ U _ 2 _ 1 _ 1  . 

1-10.24(a2-b2)+j3.2(a1-b1-10.24)1 

or 

f^(a0>a1,a2,b:j  ,b2)  =  (aQ-10. 24b2)  2+10. 24b2-104 . 8576  (a2-b2^ 2 

-  10. 24(a  -bj-10. 24)2  =  0  (6.5) 


Eqn.  (6)  is  a  set  of  high  order  simultaneous  nonlinear  algebraic  equa¬ 
tions  which  are  very  difficult  to  solve.  The  Newton-Raphson  method  that 
is  available  in  most  digital  computers  as  a  computer  program  package 


J 
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(called  the  z  systeir  [15])can  be  used  to  solve  the  nonlinear  equations. 
However,  it  is  well  known  that  the  Newton-Raphson  method  will  only  con¬ 
verge  for  a  small  range  of  starting  values  or  the  initial  guesses.  In 
order  to  improve  the  speed  of  convergence  of  the  method  four  different 
methods  of  finding  good  initial  guesses  will  be  discussed  in  the  next 


chapter. 


CHAPTER  III 

THE  INITIAL  GUESS 


In  this  report,  the  Newton-Raphson  multidimensional  method  is 
suggested  for  solving  nonlinear  equations.  However,  as  it  is  well 
known,  high  order  nonlinear  equations  have  many  solutions  and,  depend¬ 
ing  on  the  starting  values  or  the  initial  guesses,  a  solution  may  or 
may  not  be  obtained.  Therefore,  the  solution  and  the  speed  of  conver¬ 
gence  of  a  numerical  method  for  solving  nonlinear  equations  depend 
heavily  on  the  initial  guesses.  In  numerical  mathematics,  as  well  as 
in  other  areas  of  science,  finding  an  appropriate  initial  guess  for  a 
set  of  nonlinear  equations  is  itself  a  big  problem  to  be  solved.  In 
this  report,  the  following  methods  are  proposed  for  good  initial  guesses. 
The  applications  of  these  methods  depend  on  the  type  of  problem  one  lias. 


(1)  Initial  guess  by  the  synthesis  method. 

Suppose  only  the  dominant  frequency-response  data  in  (2)  are 

available  and  an  approximate  transfer  function  T  (s)  of  the  desired 

* 

Tf(s)  in  (5.1)  is  required.  The  T^Cs)  is 


*  *  *  2 
_  yyv 

Ar^s'  *  *  *  2  3 

a^+a^s+a^s  +s 


(7) 


where  a,  and  b.  are  the  initial  guesses  of  the  numerical  method.  In 
ii 

* 

the  synthesis  method  T^(s)  in  (7)  is  obtained  as  follows: 

* 

Step  1.  In  this  step  a  second-order  approximate  transfer  function  T?(s) 
is  obtained  by  using  ^  =  5.7°  and  u>c  *  3.2  rad/sec.  in  (2.6)  and  (2.7). 
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•  .38- 


* 

This  T?(s)  is 


T*  (s) 


2 

u 

n 


2  2 
s  +2c;u)  s+o) 
n  n 


G2(s) 

1+G*(s) 


(8.1) 


where  r,  =  the  damping  ratio  and  w  =  the  natural  angular  frequency. 

n 

By  following  the  basic  definitions  of  and  4>m  the  following  equations 
are  obtained. 

From  (8.1) 


G„(s)  = 


s(s+2<;u  )  2 

n  s  +2£u>  s 
n 


2;w 


G2(jw)  =  ^ 


-u)“+j2^w  co  2~  2  2 

n  /w  +4c  w  u> 


It 


180°+tan 


By  definition  jGjOw  )l  =  1,  where  wc  =  3.2  rad/sec. 


=  1 


/V,  2  2  2 

/w  +44  w  w 
c  c  n 


or. 


u4-40.9642u)2-104.8576  =  0 
n  n 


(8.2) 


Next,  by  definition 


<t>m  =  /G*(JUe)  +  L80°  =  5.7°  given 


1  2tWn 

5.7°  =  -180°  +  tan - -  +  180° 

u> 

c 


2^u> 

°r,  *  tan  5.7° 


Step  2.  In  this  step  a  third-order  transfer  function  T^Cs)  is  construct¬ 
ed  by  inserting  a  pole  (s  =  -p)  in  it  and  modifying  the  term  in  the 

•k  * 

numerator  of  T2(s)  so  that  the  steady  state  value  of  the  T^Cs)  is  equal 
to  unity,  or 


T*(s) 


2 

P  U) 

n 


2  2 
(s  +2^o)  s+to  )(s+p) 
n  n 


_ 10.290883p _ 

(s2+0 . 3194024S+10 . 290883) (s+p) 


•k 

S(s) 

l+G^s) 


(8.5) 


The  unknown  constant  p  is  determined  by  using  the  condition  in  (2.2), 

•k 

or  Re[C  (ju)  )]  =  -1.5,  where  to  =  1.9  rad/sec.  Thus,  from  (8.5) 

S  er  en 


* 

G3(s) 


_  lQ.290883p _ 

s3+(p+0.3194024)s2+(0.3194024p+10.290883)s 


(8.6) 


let  s  =  jto,  then 


\  10. 290883p 

G  (jw)  =  — . - ^ - - - - —y- 

-w  (p+0.3194024)+ju>(0.3194024p+10.290883-<o  ) 


(8.7) 


Re[G3(ju) ]  =  — 


-10.290883P  to  (p+0. 3194024) 


eT  ( p+0 . 3194024) 2+u2  (0 . 3194024p+10 . 290883-u)2 ) 2 


(8.8) 


at  to  =  io^  =  1.9  rad/sec,  RelG^tjio)]  =  -1.5 


-10 ■ 290883 (p+0. 3194024) 


3.61(p+0. 3194024)  +(0. 3194024p+6.b80883) 


-1.5  (8.9) 


After  simplification,  (8.9)  becomes 


p2  -  1.39l925823p-14. 29298735  *  0 
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or  p  =  4.540095027,  we  neglect  the  negative  root. 


Thus  (8.5)  becomes 


T*(s)  =  _ 46.72158673 _ 

3  46. 72158673+11. 74100025S+4.859497427s2+s3 


(8.10) 


Step  3.  In  this  step  another  third-order  transfer  function  (s)  is 
constructed  by  inserting  a  zero  in  (8.10)  as  shown  below. 


^  46 . 72158673+b  s  ... 

r-i  =  ~  ~  2  3  =  v 

J  46. 72158673+11. 74100025s+4. 859497427s  +sJ  1+G3  (s) 

(8.11) 

*  ... 

The  unknown  constant  can  be  determined  by  using  the  condition  in 

(6.0)  and  (2.1),  or  Re[Ge(j0)]  =  -2.1  as  shown  below.  From  (8.11),  we 


** 

C3  (s) 


b. s+46. 72158673 
**  1 

Q  (s)  =  ~  "  n  "  — — 

J  s  +4. 859497427s  +(11. 74100025-b1>s 


**  46 . 72 158673[ 1+  46-72i58673^ _ 

(s)  =  r  0^0497427 

s( 11 . 74100025-b^)  [1+  jj. 74100025-b1  s+* 


According  to  Eq.  (6.0) 


o  i  m  **  . 

Rels 


46.72158673  ,  1 _  4.859497427  . 

11.74100025-b1  '46 . 7215867  3  “  11.74100025-b^ 


Given  Re  [G  ( jO) ]  =  -2.1  in  (2.1) 
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46.72158673  bj  4.859497427 

11.74100025-b1  '46 . 72158673  "  11. 74100025-b^ 


or  b^-34.15563709b1  +  56.76713817  =  0 


which  gives  b^  =  32.4037687.  since  we  are  interested  in  the  positive  value 
only. 

Substituting  this  into  (8.11),  we  have 


T3  (s) 


46. 72 16+32. 4038s 


46 . 72 16+11 . 7410s+4 . 8595s2+s3 


(8.12) 


Equation  (8.12)  can  be  considered  as  an  approximation  of  (7)  by  assuming 
b0  =  0.  Thus  the  initial  guesses  in  (7)  are  aQ  =  46.7216,  a^  =  11.7410, 


* 

■>2 

* 


=  4.8595,  b^  =  32.4038,  and  b^  =  0.  For  solving  Eq .  (6.1)-(6.5)  these 
constants  are  used  as  initial  guesses  for  the  Newton-Raphson  method  [15). 
It  is  found  that  the  numerical  method  converges  at  the  9th  iteration 
with  the  error  tolerance  of  10  ^ .  The  solutionsof  (6.1)-(6.5)  are 
aQ  =  6.378070,  =  10.462220,  a2  =  1.259008,  b±  =  20.55667  and 

b?  =  0.243466.  Therefore,  the  desired  transfer  function  Tr(s)  is 


„  ,  N  6. 378070+20. 55667s+0. 243466s 

T  (S)  - ~2  3 

6 . 378070+10 . 462220s+l . 259008s  +s 


(9) 


The  system  represented  by  Eq .  (9)  has  the  exact  frequency  response  data 
speci f led  in  (2)  . 

(2)  Initial  guess  by  complex-curve  fitting  and  continued  fraction 

methods 

In  this  part  a  simple  method  is  presented  to  determine  the  ap- 
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proximate  coefficients  of  a  transfer  function  using  the  real  parts  and 
imaginary  parts  of  the  limited  number  of  frequency-response  data  that 
are  available.  Using  these  data  a  low-order  model  is  constructed. 

The  low-order  model  is  then  expanded  into  a  continued  fraction  of  the 
second  Cauer  form  to  obtain  a  set  of  dominant  quotients.  Some  non¬ 
dominant  quotients  are  then  inserted  into  the  continued  fraction  to 

obtain  an  amplif ied-order  model  [16],  which  is  the  desired  approximate 

* 

transfer  function  T^(s)  for  the  use  of  the  initial  guess. 

Consider  the  transfer  function 


* 

T 


r 


(s) 


b  +b  s+b.s2+. . .+b  sm 
0  1  2 _ m 

1+a, s+a0s“+. . .+a  sn 
12  n 


(10.1) 


where  a^  and  b^  are  unknown  coefficients  to  be  determined.  The  problem  of 
finding  unknown  coefficients  of  a  transfer  function  as  a  ratio  of  two 
frequency-dependent  polynomials  has  been  investigated  by  Levy  [17]. 

His  method  minimizes  the  sum  of  squares  of  the  errors  at  arbitrary  ex¬ 
perimental  points.  However,  for  finding  the  unknown  coefficients  of  a 
transfer  function  the  method  presented  next  is  comparatively  simple  and 
straightforward. 

Substituting  s  =  into  (10.1),  we  have 

,  bO+ju)kbl+(juik)2b2+"*  +  (ju,k)mbm 

Tr(J“k)  - - 2 - n - 

1+jwkal+( ju’k)  a2+. .  .+(j<«>k)  a^ 

Separating  the  real  parts  and  imaginary  parts  in  the  numerator  and  denom- 
* 

inator  of  T  (jm,  )  we  have 
r  J  k 
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2^6)  3  5  7 

,  (b0'b2a’k+b4\'b6\+-  • ' )+j  (bl\‘b3\+b5\'blV-  • 1 5 

TrU“k  =  - 2 - 4 - 6 - 3 - 5 - 7 - 

(l“a2wk+a4tVa6wk+-  *  ‘ )+j  (al‘Va3u)k+a5tVa7U)k+*  *  * ) 


=  R(u>k)  + 


‘V  +  jlk 


(10.2) 


where  and  1^  are  the  given  valuesof  the  real  and  imaginary  parts  of 

'k 

the  T^(s)  at  the  available  frequencies  Multiplying  both  sides  of 

(10.2)  by  the  common  denominator  and  separating  the  real  and  imaginary 
parts,  we  have 


(b0'b2Vb4Vb6“k+-  *  *  )+j  (bl‘Vb3U,k+b5a)k“b7Uk+'  •  *  > 


=  Rk-a2RkVa4RkVa6RkV-*-  -alWa3IkVa5Ik“k+*-* 

+  j  <aiVk'a3RkU,k+a5Rkuk"a7RkUSk+'  *  •+Ik-a2Ik\+a4IkVa6Ik\+>  •  •  > 


Equating  the  real  and  imaginary  parts  from  both  sides,  yields 


,  ,  2  ,  4,6 

b^-boW,  Tb.w.  — , co.  +. . . 

0  2  k  4  k  6  k 


Rk-a2RkWk+a4Rk\~a6RkWk+- 


alIk\+a3IkVa5Ik\+>' 


(10.3) 


bl“k‘b3“k+b5“k'b7“k+*”  =  alRku)k_a3Rkuk+a5RkVa7RkWk+- '  • 

+  Ik‘a2Ik\+a4IkVa6Iku,k+*  •  •  (10‘4) 


Eq.  (10.3)  and  (10.4)  can  be  written  as 


i 


b0'b2 Vb4\'b6\+  -  •  <+alIk\+a2iVVa3IkVa4  W*  ’  •  =  \ 

(10.5) 

3  5  7  2  3  A 

blVb3Vb5\'b7V'*’'alVk+a2Ika‘k+a3Vk-a4IkU)k---*  =  \ 

(10.6) 


In  matrix  form,  (10.5)  becomes 


1 

4 

W1 

6 

-U)^ 

Vi 

Rl“l 

T  3 
_I1U1 

“Rl“l 

i  -.22 

4 

W2 

6 

■w2  • 

I2ui2 

R2“2 

-I 

2  2 

4 

— R  a) 

2  2 

i  -u2 

4 

a*3 

6 

_a,3  * 

I3W3 

R3w3 

_T  3 

I3W3 

n  4 
-R3U3 

• 

• 

• 

• 

• 

• 

• 

• 

r 

• 

i  v 

X 

4 

w 

X 

6 

-u 

X 

I  w 

X  X 

R  J 

X  X 

<->  X 

3x 

n 

1 

-R  uj4 

X  X 

whe  re 

x  =  n 

+  1+1 

*  if 

m  is  even 

L  XJ 

(10.7) 


=  n  + 


r  if  m  is  odd 


Substituting  obtained  in  (10.7)  into  (10.6),  we  have  another  matrix 
equation  to  solve  for  b,,,  i  *  1,3,5,... 


3  5  7 

Uj  -OD^  .  .  . 


bl  1  aQ 1 1°° l+a lRlu  1 ) ” ^ a2 1 la)i+a 3*^ lu  1^ *  * '  ^ 


3  5  7  0  1  2  3 

w2  -w2  W2  -u2  **■  b3  ( ^a0I2a)2+alR2u)2^-^a2I2u,2+a3R2(‘,2^+‘  *  *  ^ 

3  5  7  . 

-w3  • '  *  °5 


3  5  7 

W  -CO  0)  —01 

_y  v  y  y 


((a_I  uP+a.K  w^)-(a„I  u)2+a.,R  u>^)+. ..) 
L_  0  y  y  1  y  y  2  y  y  3yy 
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where  =  1,  aQ=  1;  K  =  m  and  v  =  yy  if  m  is  odd;  K  =  m-1  and  y  =  -  if 
m  is  even. 

In  this  pitch  control  system,  the  available  data  is  given  in 
(2)  from  which  the  following  data  is  obtained, 


<*>!  =  wQ  -  0  ,  R2  =  TMjO)  =  1  , 


=  0 


u„  =  w  =  1.9  , 

2  en 

u  =  u  =3.2, 

3  ec 


G  (ju  )  G  (jw  ) 

r2  ■  71*2,968398’  ,3r-  >l-^»9098 

e  eir  ,e  J  en 

(ID 

G  (ju  )  G  (ju  ) 

R.  =  R e  [ T~- ]  =  0 .  3350731 ,  1-=1  [YJy-  r;-C  .-1  —  10.43159 
3  1  +G  (jw  )  3  m  1+G  (ju>  ) 

e  ec  e  ec 


Data  is  available  only  at  three  frequencies,  therefore  the  approximate 

k 

transfer  T^s)  is  assumed  to  be 


*  vbjs__ 

2  s  2 
1+n | «+n  » 


(12.1) 


Substituting  the  data  at  u^,  u2  and  u^  in  (11)  into  (10.7)  yields 


1  0 

1  -0.047786862 

1  -33.381088 


— 

—  “ 

~  — 

0 

bo 

1 

10.71591678 

al 

= 

2.968398 

3.431148544 

a2 

0.3350731 

_ 

_  _ 

__  _ 

(12.2) 


From  (12.2),  wo  get 


bQ  -  0. 047786862  a^  +  10. 71591678a2  -  2.968398  (12.3) 

b0  -  33.381088ai  +  3. 431148544a2  -  0.3350731  (12.4) 
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Substituting  bQ  =  1  into  (12.3)  and  (12. A)  yields 

-0.047786862a1  +  10. 71591678a2  =  1.968398 
-33. 381088a1  +  3. 4311A854Aa2  -  -0.6649269 

Solving  these  two  equations,  we  get 


a  =  0.0388179596 
a2  =  0.1838622891 

Then  substituting  a^  and  the  data  at  ui^  into  (10.8)  yields 

3.2b  =  9.250106342 

.\b  =  2.890658232 


Substituting  a^  and  b^,  into  (12.1)  gives 


_ 1+2. 890658232s  _ _ _ _ 

1+0. 0388179596S+0. 1838622891s 


(12.5) 


However,  the  desired  approximate  transfer  function  in  (7)  is  a  third- 

* 

order  function.  Therefore  T2(s)  in  (12.5)  needs  to  be  amplified.  In 
this  case  this  is  done  by  using  the  continued  fraction  method  [16]  as 
f ol lows . 

* 

T2(s)  is  first  expanded  into  a  continued  fraction  of  the  second 


Cauer  form  to  obtain  a  set  of  dominant  quotients.  They  are  given  as 
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h  =  1,  h2  =  -0.3506507744,  h3  *  0.9650474175  and  h4  =  16.072551656. 
Then  che  order  of  T2(s)  is  amplified  by  inserting  two  nondominant  quo¬ 
tients  hc  =  100  and  h.  =  0.1,  or 
5  o 


* 

T2(s) 


_ 1+2.890658232S _ 

1+0. 0388179596S+0. 1838622891s2 


2  1 


(12.6) 


2  1 


j*  +_L 


Substituting 


hl  “  1 

h2  =  -0.3506507744 

h3  =  -0.9650474175 

h.  =  16.07251656 
4 
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h5  =  100 
h6  =  °*1 


into  (12.6),  it  becomes 


s  .  54. 3885+162. 6914s+15. 8219s 
T.(.s)  -  T  (s)  =  T  (s;  = - r - ^ 

54. 3885+7. 5839s+10. 2146s  +s 


(12.7) 


In  solving  (6.1)-(6.5)  if  we  use  the  coefficients  in  (12.7)  as  initial 
guesses;  a^  =  54.3885,  =  7.5839,  a^  =  10.2146,  =  162.6914  and 

iV 

b2  =  15.8219,  we  have  the  desired  coefficients  in  (9)  at  the  15th  iter¬ 
ation  [15]  with  the  error  tolerance  of  10  b.  This  proves  once  again 


that  if  the  inserted  positive  quotients  >>  1  and  h^+^  <<  1  (i  is  an 
odd  number)  the  amplified  order  model  is  a  good  approximation  of  the 


original  low-order  model. 


(3)  Initial  guess  by  continued  fraction  method  [18] 


Shieh  [3]  and  Chen  [10]  have  proposed  a  continued  fraction  method  for 
model  reductions.  In  this  case  their  method  is  utilized  to  find  initial 
guesses  to  solve  Eq.  (6.1)— (6.5).  The  numerator  polynomial  N(s)  and 
the  denominator  polynomial  D(s)  in  Eq.  (1.5)  are  arranged  in  ascending 
order  and  expanded  into  the  continued  fraction  of  the  second  Cauer  form 
by  performing  repeated  long  divisions  as  follows. 


T  (s) 
e 


N(s)  b!0+b9s+b8s  +-*-+V 

"dTsT  =  7 — 7 — 27  7  IT  where  V  bi  are  glven 

all+a10S+a9S  +  ”’+V  in  (1.5) 


(13.0) 


V  V 
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where  =  1 


h?  =  -0.401749 

h3  =  -0.473321 

h,  =  25.1998 

h5  =  -0.0322195 

h£  =  -24.1061 
o 

ll  y  “  •  •  • 


(13.1) 


^22 


The  reduced  order  models  of  T  (s)  can  be  obtained  by  retaining  the  first 

e 


few  dominant  quotients,  h^  =  1,2,. 


The  number  of  quotients  used  de¬ 
pends  on  the  order  and  form  of  the  reduced  model.  This  is  explained  be¬ 
low 


1  2 

V>  1  — r  ‘  hTh’+s 

"l+  h2  1  2 


(13.2) 


V 


V 


hl+ 


V  - 


V 


h2h3+s 


hlh2h3+^hl+h3^S 


h2h3h4+(h2+h4)s 


hlh2h3h4+(hlh2+hlh4+h3h4)S+S* 


(13.3) 


(13.4) 
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h„h-h. hc+(h.h,+h„hc+h .  h_ )s+s 
2  3  A  5  2  3  2  5  4  5 


hl+ 


V 


V 


^1^2'13^4^5+^1^2^3+^1^2^5+^1^4^5+*13^4^5^  s+^l+^3+^5^ sil 

(13.5) 


V 


V 


h2+ 


h3+ 


h4+ 


h  h.h  h  h  +(h.h  h  +h,h,h,+h,hch,+h,hrh,)s 
_ 23456  234  236  256  456 _ 

h  h  h  h  hch,+(h.Lh,h,+h.h,h,h,+h,h,hth,+h,h,hch, 
123456  1234  1236  1256  1456 


+(h  +h.+h,)s 

_ 2  4  5 _ • _ 

+h.h. h.h  )s+(h  h.+h  h.+h  h,+h,h  +h,h,+hch  ) s2+s3 
3456  12  14  16  34  36  56 


(13.6) 


Substituting  the  h^'s  in  (13.1)  into  (13.6)  yields  the  third-order  ap¬ 
proximate  model  of  T  (s)  as  follows: 


3. 7 376+10. 4692S+0. 6920s2 
3. 7  376+10 . 1661s+0. 9488s2+s3 


(14) 


Using  the  coefficients  in  (14.1)  as  the  initial  guesses:  a^  =  3.7376, 
a*  =  10.1661,  a^  =  0.9488,  b*  =  19.4692  and  b*  =  0.6920,  the  desired 
solut ion  (T^ (s)  given  in  (9))of  the  set  of  nonlinear  equation  (6. 1)  —  (6.5) 
are  obtained  at  the  8th  iteration  with  the  error  tolerance  of  10 

As  it  lias  just  been  shown  in  this  particular  case,  the  continued 
fraction  method  of  finding  the  initial  guess  has  worked  out  nicely. 

However,  this  is  not  true  always.  For  example,  if  the  reduced  order  model 
by  the  continued  fraction  method  turns  out  to  be  an  unstable  system. 
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the  coefficients  of  such  a  reduced  order  model  cannot  be  used  to  solve  a 
set  of  nonlinear  equations.  Because,  an  unstable  initial  guess  often 
leads  to  solutions  which  will  give  rise  to  an  unstable  system  only.  In 
such  cases  the  following  mixed  method  can  be  used  to  obtain  a  stable 
reduced-order  model  for  approximation. 

(4)  Initial  guess  by  using  the  mixed  method. 

In  this  section  of  the  report  two  mixed  methods  are  discussed. 
One  has  the  advantages  of  both  the  continued  fraction  method  [3,10]  and 
che  dominant  pole  method  [19].  The  other  has  the  advantages  of  the 
continued  fraction  as  well  as  the  Routh  table  [11],  from  which  the 
equivalent  dominant-poles  can  be  obtained.  The  reduced-order  models 
obtained  by  the  mixed  method  are  stable  and  can  be  used  as  good  initial 
guesses . 

The  relationship  between  the  quotients  h^  and  the  coefficients 
and  b^  in  (13.0)  can  be  expressed  by  the  following  matrix  Eq.  [3,4]: 

[b]  =  [H]  [a]  (13) 

T 

where  [a]  =  [  a  j  ,^n_2 >  •  •  •  .  a2  ,a]L  ,aQ]  , 

=  [bn-l,bn-2’",,b2’bl,b0]  * 

[HJ  =  [H2 ]_1  [Hjl  , 


here  T  designates  transpose  of  a  matrix 
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I 

i 


IHJ  = 


h  0  0 


1  h2  0 


0  0 


0  0  0 


[H, 


0  0 


o  h2  0 


0  1 


0  0 


0  0  0 


0  0 
0  0 
0  0 
0  0 


1  h 


1  0 
o  h. 


0 

0 


0  1 


0  0 


0  0  0 


0  0 
0  0 
0  0 
0  0 


1  h 


10  0 
0  1  0 

0  0  h„ 


0  0 


0  0  0 


0  0 
0  0 
0  0 
0  0 


1  h 


n-i 


10  0 
0  1  0 
0  0  1 
0  0  0 


0  0  0 


0  0 
0  0 
0  0 
0  0 


1  h 


n-1 


10  0 
0  1  0 
0  0  1 
0  0  0 


0  0  0 


Consider  the  reduced-order  model  of  the  original  system  as 


T  (s)  = 
r 


e_+e,s+...+e  ,s 
0  1  r-1 


r-1 


d_+d  s+.,.+d  ,sr  1+d  sr 
0  1  r-1  r 


0  0 
0  0 
0  0 
0  0 


0  h 


0  0 
0  0 
0  0 
0  0 


0  h 


(16) 


The  denominator  polynomial  in  (10)  is  approximated  by  the  pro¬ 


duct  of  the  dominant  poles  of  the  original  system  Tg(s).  Thus  is 


known.  Replacing  and  in  (15)  by  d^  and  e^  in  (16),  respectively, 


Eq.  (15)  can  be  solved  for  e ^  in  (16).  The  T^(s)  obtained  has  the 


dominant  poles  and  the  dominant  quotients  of  ^(s)  and  it  is  always 


stable,  therefore,  T  (s)  can  be  used  as  a  good  initial  guess  in 


solving  (6.1)-(6.5). 

In  case  the  roots  of  D(s)  in  (13.0)  are  not  available,  the  ap- 
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proximate  equivalent  dominant  poles  and  the  resulting  coefficients  d. 
can  be  determined  from  the  Routh  table  as  suggested  by  Hutton  and 
Friedland [ 9 ]. The  steps  involved  are  explained  below. 


Assume  T(s) 


,  n— 1  ,  n-2  , .  .  , 

b  ..  s  +b  „s  +. . .+b„+b,+b„ 
n-1 _ n-  2 _ 2  10 

n_l_  n_l,  ,  . 

a  s  +a  ,  s  +. . .+a„+a,+a_ 

n  n-1  210 


n  ( s) 
d(s) 


(17.1) 


is  the  original  transfer  function  for  which  the  reduced  order  model  is 
needed . 


Step  1.  Construct  a  Routh  array  [20]  using  the  coefficients  a^  of  the 
d(s)  above  and  the  Routh  algorithm.  The  Routh  array  is  shown  below. 

To  obtain  a  general  algorithm  a^  is  expressed  double-subscripted  nota¬ 
tion,  for  example,  a.  . 

i.l 


A  A  A 

—  ci  3.  a  s  a  r)  a.  ~  -  a  . 

n  12  n-2  13  n-4 


~  a  i  =  3  q  3a a  ~  a  c 

n-1  22  n-3  23  n-5 


a12_Yla22  a32  '  a13"Yla22  a33 


a22-Y2a32  a42  =  a23_Y2a33 


a 


0 
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I 

I 

I 


n-2,2 


Vl,2  =  a0 


=  a. 


In  general 


ai,j  "  ai-2,j+rYi-2ai-i,j+i;  1  1,2 . j  3>4, 

'i  =  3i,l/ai+l,l 


Stej£  J2 .  In  this  step  various  approximate  low-order  polynomials 
constructed  from  any  two  consecutive  rows  in  the  Routh  array, 
say  from  the  last  row  and  the  next  to  the  last  now  and  so  on. 
explained  below. 

The  first  order  (i  =  1)  approximate  equation  is 


* 

d^s) 


a  ,  s+a  =  a  s+aA  =  0 

n,l  n+1,1  n,l  0 


The  second  order  (i  =  2)  approximate  equation  is 


*  2  2 
d„(s)  =  a  .  ,s  +a  .s+a  -  a  .  ,s  +a  ,s+a„  =  0 

2  n-1,1  n,l  n-1,2  n-1,1  n,l  0 


(17.2) 


(17.3) 

it 

d^ts)  are 
for  example, 
This  is 


(17.4) 


(17.5) 
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The  third  order  (i  =  3)  approximate  equation  is 


*  3  2 

d..(s)  =  a  -  .s  +a  s  +a  .  „s+a  .  =  0 

3  n-2,1  n-1,1  n-2,2  n-1,2 


3  2 

=  a  „  ,s  +a  .  ,s  +a  „  .s+a.  =  0 
n-2,1  n-1,1  n-2,2  0 


(17.6) 


and  so  on. 

When  the  original  system  (17.1)  is  asymptotically  stable,  all  y^are  posi- 

k 

tive  values  and  the  approximate  polynomials  d^(s)  are  the  Hurwitz  poly- 

■k 

nomials.  The  d.(s)  are  normalized  simply  by  dividing  each  coefficient 

*  .  j  * 

in  d^(s)  by  the  coefficient  of  the  highest  order  term  in  s.  These  normalized  d^(s 

•k 

are  the  denominator  polynomials  of  the  reduced-order  models  T^(s)  of  the 

* 

original  system.  Then  the  numerator  polynomials  of  T.(s)  are  determined 

* 

simply  by  substituting  the  coefficients  of  d^(s)  in  place  of  [a]  in  (15) 

norm 

and  then  solving  the  matrix  equation  (15)  for  [  b] ,  which  are  the  coeffi- 

k 

cients  of  the  numerator  polynomial  of  T  (s) .  * 

1  *  n3(s) 

The  third-order  reduced  order  model  T.  (s)  =  — r -  of  the 

3m  j  e  \ 

d3(s) 

original  pitch  control  system  in  (1.5)  obtained  by  using  the 

mixed  method  is  explained  below. 

At  first,  the  Routh  array  of  the  pitch  control  system  in  (1.5) 

is  obtained.  From  the  Routh  array  the  normalized  approximate  denominator 
* 

d3(s)  is  found. 


d*(s)  =  sJ+0.9523822967s2+10.19241445s+3. 745517989 


(17.7) 


*  A  3  * 

To  determine  n3(s)  =  b2S  +b^s+bg,  the  coefficients  of  d3(s)  are 


substituted  into  (15)  as  shown  below. 
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where  the  h^'s  are  the  quotients  of  Tg(s)  in  (1.5),  which  are  given  in 
(13.1).  Substituting  the  values  of  h^,  and  h^  from  (13.1)  into  (17.9) 
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and  then  simplifying,  we  get 


bQ  =  3.7455  (18.1) 

bfl  -  0.401749  b  =  -4.094787 

substituting  (18.1)  yields 

b  =  19.5154  .  (18.2) 

and  0.524679  b  +  0.19096  b2  =  10.37427 

Substituting  (18.2)  yields 


b9  =  0.7066 


(18.3) 


Therefore 


0.7066s  +19. 5154s+3. 7455 


s3+0.9524s2+10.1924s+3.7455 


(19) 


In  solving  the  nonlinear  Eqs.  (6.1)-(6.5)  if  the  coefficient  of  T^^s) 


*  * 

in  (19)  are  used  as  starting  values:  =  3.7455,  a^  =  10.1924, 


a2  =  0.9524,  b^  =  19.5154  and  b2  =  0.7066;  the  Newton-Raphson  method 
[151  converges  to  the  desired  solution  in  (9)  or 


t  (  \  6. 37807+20. 55661S+0. 243466s2  Gr(s) 

Ir's'  2  3  “  1+G  (s') 

6. 37807+10. 46222s+l. 259008s  +s  rv  ' 


at  the  8th  iteration  with  the  error  tolerance  of  10 


1 
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From  (9) 


G 

r 


(s)  =  The  open-loop  transfer  function  of  the  standard  model 


T 

r 


(s)  . 


6. 37807+20. 55661s+0. 24346s2 
s(-10. 09439+1. 015542s+s2) 


(20) 


The  Nyquist  plot  of  Gr(s)  is  shown  in  Fig.  2  and  the  unit  step 
responses  of  T^Cs),  T3((s),  T3m(s)  and  Tg(s)  are  compared  in’Fig.  3.  All 

"ft  if 

three  reduced-order  models  T^Cs),  T-^Cs)  and  T^s)  give  very  satisfactory 
approximate  time  response  curves.  However,  only  the  T^(s)  in  (9),  which 
uses  the  method  of  dominant  frequency  response  data  matching,  has  the 
exact  dominant-frequency  response  data  as  the  original  system  T^(s) 
given  in  (2). 


ORIGINAL  11TH  ORDER  SYSTEM:  T' 


™E  RE$P0NSES  °F  *>  ™.R0  ORDER  REDUCED  «0EL5. 


CHAPTER  IV 


SIMPLIFICATION  OF  THE  EXISTING  STABILIZATION  FILTER 


As  it  appears  from  its  name,  the  purpose  of  the  stabilization 

filter  is  to  stabilize  the  original  unstable  system.  The  transfer 

function  of  the  existing  stabilization  filter  F  .  (s)  is  known  and  is 

stflb 

given  in  (1.2).  As  it  is  mentioned  in  the  introduction  of  this  report, 
the  objective  of  this  report  is  to  redesign  the  stabilization  filter 
so  that  the  cost  of  implementation  can  be  reduced  and  at  the  same  time 
the  performance  of  the  redesigned  pitch  control  system  is  the  same  as 
that  of  the  existing  stabilized  pitch  control  system. 

In  this  chapter  two  different  transfer  functions  are  obtained 
for  the  stabilization  filter.  Both  of  these  transfer  functions  are  ob¬ 
tained  by  direct  simplification  of  the  available  transfer  function  of 
Fstab^S^’  and  °ne  t'lera  *s  obtained  by  using  the  dominant  data  match¬ 
ing  method  of  Chapter  II. 

The  F  ,  (s)  in  (1.2)  can  be  considered  as  the  closed-loop 
stab 

transfer  function  of  a  control  system  as 


N  (s)  G  (s) 

F  /  v  _s_ _ _  _ stab  _ 

stab  S  D  (s)  1+G  ,  ( s ) 

S  SlaO 


460800s2+69120000s+144x107 

s4+250s3+76900s2+72x105s+9*108 


(21.1) 


where  the  open-loop  transfer  function  G  ,  (s)  is 

stab 


G  (s)  _ 460800s 2+6 9 12QQ00s+ 144* IQ7 _ 

stab  s4+250s3-383900s2-61920000s-5.4*108 


(21.2) 


The  dominant  frequency-response  data  of  this  system  are  given  below. 
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i) 

Gstab(j0)  =  "  0.375 

(22.1) 

ii) 

Re[Gstab(^sl,>l  -  -1.032833 

(22.2) 

Im[G  „  .  (ioj  )]  =  0.002017351 
stab  sw 

(22.3) 

where 

at  =  The  phase  crossover  frequency  of 

=  140  rad/sec. 

the  stabilization  filter 

iii) 

R‘IC,t.b«».c)1  ■  -1-002941 

(22.4) 

I",ICstab<1"sc)1  "  -°-<>“68759 

(22.5) 

where  wgc  =  The  gain  crossover  frequency  of  the  stabilization  filter 
=  200  rad/sec. 


Suppose  the  reduced-order  model  F  ,  (s)  of  the  stabilization  fil- 

sl 


ter  is 


F  (s)  = 
si 


w 


Gsl(s) 


a^+a^s+s 


(23.1) 


where 


Csl<s> 


The  open-loop  transfer  function  of  Fg^(s) 


vv 


(a0-b0)+(a1-b1>s+s 


(23.2) 


The  constants  a^  and  b^  are  unknown  constants  to  be  determined.  Using 


the  specifications  given  in  (22)  and  following  the  basic  definitions  of 


those  specifications  the  unknown  constants  and  b^  are  determined  as 
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shown  below. 

For  F  ,  (s)  in  (23.1)  to  be  a  reduced  order  model  of  F  K(s) , 
sl  staD 

G  „(s)  must  satisfy  all  the  specifications  of  G  ,  (s)  in  (22).  Apply- 
sl  stab 

ing  the  condition  in  (22.1)  to  the  system  Gg^(s)  in  (23.2)  yields 


at  s  =  jO  Gg^(jO) 


or, 


Vbo  °’375 


b0  =  ll6ao 


(24.1) 


Substituting  (24.1)  into  (23.1)  and  (23.2),  respectively,  we  get 


and 


rsl(s)  ■ 


1. 6a  +b  s 
2 

a^+a^s+s 


(24.2) 


Gsi(s)  = 


1 . 6a^+b^s 


-0.6a0+(a1~b^)s+s 


(24.3) 


at  s  =  jui  Ggl(ju)) 


1 . 6aQ+jtub^ 


-(0.6aQ+a)  )+jw(a1~b1) 


(1.6aQ+jub1)[-(0.6aQ+w  J-juKa^-bj)  ] 
(0. 6aQ+o)3)7+(D7  (a^-b^)^ 


-1.6an(0.6an+u7)+o>^b  (a..-b  ) 

.-.  Re[G  .(ju)]  =  - - - 5-5—5 - - 

(0.6an+w  )  +10  (a^-b^) 


(24.4) 


0 


Im^iO)  1 


-a)bj  (0. 6aQ+w  )-l .  6waQ (a^-b^) 
(0. 6aQ+u>7)  (a^-b^) 


(24.5) 


i) 


Specification  in  (22,2)  yields 


44 


Re[G  (jl40)]  =  -1.032833 

Substituting  (24.4)  above  gives 

-1.6a  (0.6a  +19600)+19600b  (a  -b  ) 

- y - y — - — L_  =  -1.032833 

(0.6aQ+19600)  +19600(a1-b1) 

or  f^taQ.a^.b^)  =  -1.6a0(0.6a0+19600)+19600b1(a1~b1) 

+  1. 032833 [  (0.6a0+19600)2+19600(a1-b1).2]  =  0 

(25.1) 

ii)  The  data  in  (22.3)  when  applied  to  (24.5)  yields 

Im[Gsl(jl40)]  =  0.002017351 

-140b  (0.6a  +19600)-224a  (a  -b  ) 

or,  - i - - - y — 11_  ,  0.002017351 

(0. 6a0+19600)z+19600(a1-b1) 

or  f2(a(),a1,b1)  =  -140b1(0.6a0+19600)-224a()(a1-b1) 

-  0. 002017351[(0. 6aQ+19600)2+19600(a1-b1)2]  =  0 

(25.2) 

iii)  The  data  in  (22.4)  when  applied  to  (24.4)  gives 

Re[Gsl(j200)l  =  -1.002941 

-1 . 6a  (0. 6an+40000)+40000b.  (a^b,  ) 

_ o _ y _  iii 

(0. 6a0+40000) 2+40000(a1-b1) 2 


or 


-1.002941 
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or  f ^(a^.a^ ,b^)  =  - 1 . Sa^ (0. 6aQ+40000)+40000b^ (a^-b^) 

+  1.002941[ (0.6a0+40000)2+40000(a1-b1)2]  =  0 


(25.3) 


Equation  (25)  is  a  set  of  nonlinear  equations.  The  unknown  con¬ 
stants  ai  and  b^  in  (23.1)  are  determined  by  solving  (25).  However,  to 
solve  equation  (25)  the  proper  initial  guesses  have  to  be  determined  first. 
As  discussed  in  Chapter  III,  the  initial  guesses  can  be  determined  from 

the  reduced-order  model  of  the  existing  stabilization  filter  F  .  (s) 

stab 

in  (1.2).  Using  the  mixed  method  of  the  continued  fraction  approximation 

* 

and  the  Routh  approximation,  a  reduced-order  model  Fr^(s)  is  obtained  as 


follows.  *  * 

*  b„+b,s  *,  s 

Assume  F^(s)  =  j  = 

aQ+a^s+s  d  (s) 

Routh' s  st  ilitv  criterion  for  F  _  ,(s)  is 

stab 


(26.1) 


i  48100 


2522245.322  ' 

-  —  -  -  _  _  j 


76900 


72x10' 


As  discussed  in  the  Section  4  of  Chapter  III, the  d  (s)  in  (26.1) 
is  approximated  from  the  Routh  criterion  shown  above.  Thus 


d*(s)  =  48100s2+2522245.322s+9*108  *  0 


After  normalization  d  (s)  becomes 


d  (s)  =  s  +52. 4375s+18711. 01871 


Therefore,  now  (26.1)  becomes 


Frl(s)  =  -j 


•k  -k 

bis+bo 


s  +52.4 375s+18711. 01871 


The  quotients  h.  of  F  ,  (s)  are  obtained  below 
l  stab 


72x10  76900  250 


144x10  69120000  460800 


-36x10  -211100  250 


0676000  470800  40 


?V  k 

Using  Eq.  (15)  b^  and  b^  in  (26.1)  are  determined  as  follows 


~hl 

0 

~  *  "1 

bo 

i 

0 

18711.01871 

1 

hlh2_ 

Lb;j 

0 

52.84375 

or 

— * 

0.625 

0 

r  *  - 
b0 

1 

0 

18711.01871 

* 

S 

1 

-25 

b. 

0 

-40 

52.4375 
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,*  18711.01871  - -  - - 

or  b0  oT623 -  299  37.62994 


(26.4) 


and  b*  -  25b*  «  -2097.5 


or 


*  =  2097-5  b0  =  2097.3+29937.62994 

bl  -25  25 


or  b  =  1281.40525 


(26.5) 


Substituting  (26.4)  and  (26.5)  into  (26.3)  yields 


Frl(s) 


1281. 40525s+29937. 62994 
s2+52.4375s+18711. 01871 


(26.6) 


Using  the  coefficients  of  F  (s)  in  (26.6):  a^  =  18711.01871, 

V;  * 

a^  =  52.2,V7S,  =  29937.62994  as  initial  guesses,  the  nonlinear  equa¬ 

tions  in  (25)  are  solved  by  the  Newton-Raphson  method  [15]  and  the  fol¬ 
lowing  solutions  are  obtained  at  the  7th  iteration  with  the  error 
tolerance  of  10  b : 


aQ  =  20917.459536 
a  =  29.981293 
b  =  957.260014 

Since  b  =  1 . as  in  (24.1)  b^  becomes 
bQ  =  33467.93525 

F  , (s)  the  desired  low-order  stabilization  filter  in  (23.1)  is 
s  i 
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957. 260014S+33467. 93525 
s2+29.981293s+2091 7. 459536 


(27) 


The  unit  step  response  of  the  existing  stabilized  pitch  control  system 
in  Eq.  (1.5)  and  the  redesigned  pitch  control  system  using  F^(s)  in 
(27)  and  the  Gq(s)  in  (1.3)  are  shown  in  Fig.  4.  The  result  is  fairly 
satisfactory. 

An  alternate  approach  for  redesigning  the  stabilization  filter 
by  direct  simplification  of  the  existing  stabilization  filter  is  pro¬ 
posed  as  follows: 

As  it  is  mentioned  at  the  beginning  of  this  chapter, the  function 

of  the  stabilization  filter  is  to  convert  the  dominant  data  at  u  =  0, 

w  =  1.9  rad/sec.  and  w  =  3.2  rad/sec.  of  the  original  unstable  svs- 
en  ec 

tem  Gq(s)  in  (3)  to  the  assigned  dominant  data  of  G&(s)  in  (2).  Taking 
advantage  of  this  fact,  we  can  directly  apply  the  dominant-data  matching 
method  to  fit  a  low-order  stabilization  filter  that  satisfies  the  speci¬ 
fications  assigned  in  Eqn.  (4).  Let  us  assume  that  the  desired  low-order 

model  of  F  ,(s)is 
stab 


fs2(s) 


yy 

a^+a^s+s^ 


(28.1) 


Applying  the  condition  in  (4.1)  to  F  „(s)  in  (28.1)  yields 

sz 


Re[F  _( JO) ]  -  —  =  1.6 

sZ  aQ 


b0  = 


(28.2) 


Substituting  (28.2)  into  (28.1)  gives 


THE  EXISTING  SYSTEM  Te(s)  AS  VEIL 
REDESIGNED  SYSTEM  USING  Fs?(s) 


Figure  4.  Time  Responses  of  Various  Models 
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1.6a  +bs 

5,2 <s>  '  — — 77 

a^+a^s+s 


(28.3) 


At  s  =  ju) 


1.6a  +jub 

F.,(j«)  =  - J - ~ 

(aQ-u)  )+ji»)a^ 


lFs2(jw)l  = 


Jl 


2  2  2 
2.56a‘+u 


/  2.2  2  2 

/  (aQ~W  )  +W  3 


(28.4) 


and 


/Fs2(ja))  =  tan  176^ 


-1  Wbl  -1  W31 

-  tan  - ; 


a„-w 

0 


=  tan 


_j  «b1(a0~u)  ) -1 . 6toa^a^ 
2  2 

1.6a0(a0~u)  )+w  a.jbj 


(28.5) 


At  s  s  1m  =  il.9  the  values  of  IF  0  ( j  co)  |  and  /F  ~(ju)  in  (28.4)  and 
•  en  1  s2  /  s/ 

(28.5)  respectively  are  matched  to  the  corresponding  values  of  | Fstab  ( j  1  •  9)  |  to¬ 
gether  and  /f  . (il.9)  in  (4.3).  Thus,  we  have 
j/^stab 


>/2.5^+3-61b5 

|F  2(jl.9)|  «  — — - - -  =  1.605107127 


/(a0-3.61)2+3.61a^ 


or 


f1(a0,a1,b1)  =  2.56a2+3.61b2-2.576368889[(a0-3.61)2 


+  3.61a*]  =  0 


(29.  1) 


,  1.9b, (an-3.61)-3.04ana  0 

-  1.6,1„(,!-3.61H3.61^7  - 


*0'“0 


1“1 
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or 


f0(aQ,a1,b1)  =  l^b^a^-bD-S^a^ 


0. 0759963811 [1.6a()(a0-3.61)+3.6la1b1]=0  (29.2) 


When  s  =  jw  =  j  3 . 2  the  val 


ilue  of  /. 

ec  /  s2  ec 


F  (jw  )  in  (28.5)  is  compared  with 


the  value  of  /F  ,(iw  )  in  (4.5).  Thus,  we  have 
/  stab  J  ec 

3.2b  (a  -10.24)-5.12a  a 

/Fs2(j3‘2)  =  tan  l76T0(a0-10.24)+10.24a1b1  “  7* 


293349493° 


f3(a0,a1,b1)  =  3.2b1(aQ-10.24)-5. 12aQa1 


0.1279849782[1.6a  (a  -10.24)+10.24a  b  1=0  (29. 3) 


Using  the  initial  guesses  obtained  in  (26.6)  the  set  of  nonlinear  equa¬ 
tions  in  (29)  is  solved  for  the  unknowns  a^,  a^  and  b^  by  using  the  Newton- 
Knpbson  method.  The  solutions  are  obtained  at  the  9th  iteration  with 
the  error  tolerance  of  10  The  solutions  obtained  are 


aQ  =  13301.999297 
a:  =  3.318051 
b}  =  856.628596 


Thus,  the  desired  low-order  model  in  (28.3)  is 


,  856. 628596s+21283. 19886 

t  „(s)  =  —= - 

s  +3. 318051s+13301. 999297 


(30) 


The  unit-step  response  of  the  existing  stabilized  pitch  control 


system  T£(s)  in  (1.5)  and  the  redesigned  system  that  uses  the  low-order 
filter  in  (30)  and  Gq(s)  in  (1.3)  are  shown  in  Fig.  4.  The  re¬ 

sult  is  perfect.  Comparing  the  unit-step  response  curves  in  Fig.  4,  it 
is  clear  that  as  far  as  the  performance  of  the  entire  pitch  control  sys¬ 
tem  is  concerned  F  ~(s)  in  (30)  is  a  better  filter  than  F  .  (s)  in  (27). 
s2  si 

This  implies  that  the  existing  stabilization  filter  F  ,  (s)  in  (1.2) 

stab 

might  be  overdesigned.  Obviously,  the  implementation  cost  of  the  filter 
F^Cs)  is  *ess  t*lan  that  of  Fgtab(s)  in  (1.2). 


CHAPTER  V 


REDESIGN  OF  THE  STABILIZATION  FILTER  BY  AN  ALGEBRAIC  METHOD 

In  Chapter  IV  of  this  report  the  original  fourth-order  stabi¬ 
lization  filter  Fsta^^s)  has  been  simplified  to  two  second-order 
filters,  Fg^(s)  and  F^Cs)*  using  the  dominant-data  matching  method 
discussed  in  Chapter  II.  It  is  noticed  that  all  three  stabilization 
filters,  the  original  as  well  as  its  simplified  models  consist  of  com¬ 
plex  poles.  It  is  also  observed  that  all  three  filters  mentioned  above 
are  placed  in  the  feed  forward  loop  and  as  a  result  the  system  becomes 
very  sensitive  to  external  disturbances.  If  alternate  filters  can  be 
designed  and  placed  in  both  feed  forward  and  feedback  loops, i)  the  de¬ 
signed  filters  may  turn  out  to  be  simple  transfer  functions  with  posi¬ 
tive  real  roots  and  because  of  this  it  may  be  possible  to  synthesize 
the  filters  using  passive  elements,  and  ii)  the  performances  of  the 
designed  system  can  be  greatly  improved.  The  fact  that  the  fixed 
configuration  of  the  compensators  in  the  feedback  loop  enables  the  de¬ 
signed  system  to  be  insensitive  to  the  parameter  variations  and  modeling 
errors  will  reduce  the  effects  of  external  disturbances  and  improve 
the  stability  of  the  system.  Thus  the  redesigned  feedback  system  has 
all  the  advantages  [14]  of  feedback  control  systems. 

In  this  chapter  the  algebraic  method  proposed  by  Shieli  [3]  and 
Chen  [/( ]  is  extended  and  modified  to  redesign  the  pitch  control  system 
The  steps  involved  are  summarized  as  follows. 

Step  1.  Assign  the  design  goals  using  frequency-domain  specifications 
and  model  a  standard  transfer  function,  known  as  the  standard  model, 
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using  die  dominant  data  matching  method  discussed  in  Chapter  11  of  this 
thesis . 

Step  2.  Expand  the  standard  model  obtained  in  Step  1  into  a  standard 
fraction  expansion  of  the  second  Cauer  form  by  performing  repeated  long 
divisions  as  shown  in  (13.0)  to  obtain  the  dominant  quotients.  Using 
these  quotients  obtain  the  matrix  l H ]  in  Equation  (15). 

Step  3.  Assume  the  fixed  configuration  of  compensators  with  unknown 
parameters  and  determine  the  overall  transfer  function  of  the  system. 

Thus,  the  overall  transfer  function  of  the  system  will  consist  of  the 
unknown  parameters. 

Step  4.  Substitute  the  coefficients  of  the  overall  transfer  function  ob¬ 
tained  in  Step  3  into  the  vectors  [a]  and  [b]  in  Equation  (15)  and  ex¬ 
pand  tin-  matrix  equation  (15)  to  obtain  a  set  of  equations. 

Step  5.  Solve  the  set  of  equations  obtained  in  Step  4  to  determine  the 
unknown  constants  assigned  in  the  compensators. 

The  designed  system  obtained  by  using  the  algebraic  method  has 
the  exact  dominant  quotients  of  the  standard  model.  In  other  words, 
the  designed  system  is  a  good  approximation  of  the  standard  model  that 
has  the  exact  assigned  dominant  data. 

It  is  noticed  that  the  original  unstable  system  Gq(s)  in  (1.3) 
is  a  high  order  transfer  function  with  large  coefficients.  Therefore, 
in  order  to  simplify  the  procedure,  before  proceeding  to  design  the 
pitch  control  system  by  using  the  algebraic  method,  a  reduced-order 
mode)  of  Gq(s)  is  determined  by  using  the  dominant~data  matching  method- 

1 

i 
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The  unstable  transfer  function  Gg(s)  in  (1.3)  can  be  decomposed 
into  a  stable  function  and  an  unstable  portion  as  follows: 


G0(s)  =  s (s-2 .921)  T0(S) 


(31.1) 


where  the  stable  portion  Tq(s)  is 


,  ,  324  332. 316 (s+0. 1933)  (s+65)  (s-t-15Q0) 

0(S)  (s+3. 175) (s+87 .9±j95 .5) (s+112 .5) (s+1385) 


(31.2) 


The  pole  at  the  origin  and  the  unstable  pole  at  s  *  2.921  are  consider¬ 
ed  to  be  the  dominant  poles  of  the  system.  Therefore,  they  are  retained 

•k 

in  the  simplified  model  Gq(s)  of  GQ(s) ,  or 


G0(s)  =  G0(s)  ”  s(s-2.921)  "0 


T*(s) 


(31.3) 


Whore  Tq(s)  is  the  reduced-order  model  of  Tq(s)  obtained  by  us¬ 
ing  the  dominant-data  matching  method.  The  frequency  response  data  of 
Tq(s)  that  are  used  as  dominant  data  for  the  transfer  function  fitting 

aregain  margin,  phase  margin,  phase-crossover  frequency,  gain-crossover 

k 

frequency,  and  the  final  value  at  w  =  0.  The  Tq(s)  obtained  is 


N  496.854897s  +192897. 961011s+37103. 33375 
T0(s)  =  — - 5 - 


(31.4) 


s  +117. 073733s  +16552. 300003s+50595. 685093 


The  Ty(s)  obtained  is  a  low-order  model  of  Tq(s)  with  smaller 
coefficients.  Thus,  the  design  process  can  be  greatly  simplified. 


Therefore*;^)  =  —  - 


496.854897s  +192897. 961011s+37103. 33375 


sJ+114. 1527 33s4+16210.32763s3+2246.41679s2-147789. 9961s 

(31.5) 


Following  the  steps  proposed  at  the  beginning  of  this  chapter 


the  first  step  to  design  a  system  by  the  algebraic  method  is  to  determine 

the  standard  model.  In  this  case  the  standard  model  T  (s)  has  been  de- 

r 

termined  earlier  in  Chapter  III  and  is  given  in  (9).  Writing  T^s)  once 
again  and  expanding  it  in  a  continued  fraction  expansion  yields 


6.3780  7+20 . 55 66 ls+0. 243466  s2 
6.3780 7+10.46222s+1.259008s2+s3 


where 

hl  "  1 

h2  =  -0.631845015 
h3  =  -0.476189214 

h.  =  14.799589050  (31.6) 

4 

h5  =  -0.102867450 
hfe  =  -13.924278040 

In  the  next  step  a  series  compensator  G^(s)  and  a  parallel  compensator 
G2(s)  are  assigned  as  shown  in  the  block  diagram  of  Fig.  5-1.  G^(s) 


Figure  5.  The  Block  Diagrams  of  the  Redesigned  System  Using 
Algebraic  Method 


1 
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and  G,,( 


and 


in  Fig. 


where 


s)  are  assumed  with  unknown  parameters  x^ ,  i 


1,2,. . . ,7  as 


G1(s) 


X  s+x, 
6 _ (_ 

s+xc 


(32.1) 


x.s  +x, s+x 

<*•>  ■  4-t-1 

s  +x^s+x2 


(32.2) 


The  overall  transfer  function  T^(s)  of  the  feedback  system  shown 
5-1  is 


Tf  (s)  = 


t>Q+b^s+  ...  +b^s 

- 

ag+a^s+  . . .  +3gS 


(32.3) 


aQ  =  37103.  33375x2x-, 

ax  =  192897. 961011x2x7+37103.33375(x2x6+x4x7)-]A7789.996lx2x5 

a„  =  496 . 85A897x_x_+192897 . 961011 (x_x+x  x., ) 
l  2  /  2  b  /*  / 

+  37103.33375(x.  x,+x.,x.,)+2246.  41679x  xt 
4  6  3  7  2  5 

-  147789. 9961(x2+x1x5) 

a„  =  496 . 854897 (x0x,+x, x.J+192897 . 961011 (x. x.+x,xj 

3  2647  463/ 

+  37103. 33375x_x.-147789 .9961 (x,+xc) 

3  6  ]  5 

+  2246. 41679 (x2+x^x^)+162 10. 32763x2x^ 

a,  =  496. 85489 7(x. x,+x_x.,)+192897 . 96 101 lx, x, -147789 .9961 

4  4  6  3  7  3  6 

+  2...6.41679(x]+x5)+16210.32763(x2+x]x5)+114.152733x2x5 

a,  =  496. 854897x,x,+2246 .4167 9+16210. 32 763 (x.+xc)  (32.4) 

5  3  6  1  5 

+  114. 152733(x2+x^x^)+x2x^ 

a,  =  16210. 32763+114 .152733(x,+x. )+x„+x, x. 

6  1  5  2  i  0 
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a?  =  114.152733+x]+x5 
a8  =  1 

bQ  =  37103. 33375x2x? 

b.  =  192 89 7. 96 101  lx „x ..+  37103.  33375 (x»x£+x  x_ ) 

1  Z  f  L  o  1  / 

h„  =  496 . 854897x.x_+192897 . 961011 (x„x,+x x. ) 
z  l  /  Zb  i  / 

+  37103.  33375(x.,x,+x-,) 

lb  7 

b_  =  496.854897 (x_x,+x  x7 )+l 9289 7. 961011 (x  x  +x7)+37l03. 33375x. 

3  2617  167  6 

b.  =  496.854897 (x  x,+x  )+192897 . 96 101 lx, 

4  1  o  7  6 

b,  =  496 . 854897x 
-)  o 

b6  =  o 


In  order  to  match  the  seven  unknown  parameters,  x^,  i  =  1,2,..., 7 
in  (32)  for  this  type  '1'  system  we  need  eight  quotients  h^,  i=l,2,...,8 
in  (9).  Therefore,  the  third  order  standard  model  in  (9)  with  the  quo¬ 
tients  h.  given  in  (31.6)  has  to  be  amplified  to  a  fourth-order  system. 
This  is  done  by  inserting  h7  ^  100  and  hQ  =  0.1  as  shown  below. 

,  s  6. 37807+20. 55661s+0. 243466s2 

T  (s)  = - 2  3 

6. 37807+10. 46222s+l. 259008s  +s 
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1 


=  T  (s)  = 
a 


63.780980O7+211.8989926s+22.87561717s2+Q. 34346s3 
63. 78098007+110. 9545223s+23.00917551s2+11. 30110515  s3+s* 


It  [15]  has  been  shown  that  T  (s)  in  (33)  is  a  good  approximation 

3 

of  the  original  standard  model  T^fs)  in  (9). 

Substituting  the  h^,  i  =  1,...,6  in  (31.6)  including  h^  =  100  and 

hu  =  0.1,  the1  matrices  [11,]  and  [11-  |  in  (5)  are  obtained  next, 
o  1  L 


11.301105  23.009176  110.95452  63.78098 
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Substituting  the  unknown  constants  a^  ,  i  =  0,1,... ,7,  b^,  i  =  0,1,. ..,7 
and  IH^]  and  [H0]  obtained  above  into  (15)  yields  a  set  of  equations  in 
terns  of  a.  and  b.  as  follows. 

t  l 


lb  1  =  l H  j la] 

=  ( H2 1~ 1  [Ht]  [a] 

or  [H2l[b]  =  [Hjla] 


or  expanding  the  above  matrix  equation  yields 


f 


1 


b.) 

i 


0  0 


0 


f2(VV  =  ^q-^*  6318422396  (1^-3^)  =  0 

f „(a  .  ,b  , )  =  a  +0 . 30087727  (a  -b.)-0. 52380951b  =  0 

dill  22  1 

f.(a.,b.)  =  b,+7.1203141b.+4.4528546(b.-a.1)-14. 167729a.  =  0 
4  i  i  1  2  J3  2 

f5(ai,bi)  =  a2-1-1565286a3-0.45805621(a4-b4)-0.42094152b2 
+  0.4315751b3  =  0 

f , (a . ,b  . )  =  b.+l. 2 590 lib .+10. 462223b. +6.  37 8098 (b. -a. ) 
6ii2  3  4  55 

-  0. 2 4 346000a _-20. 55667a,  =  0 
3  4 


f  (a , ,b  .  ) 
7  i  i 


a_+23.  189471a  ,+2055. 2089a,. +637 . 8098(a,-b, ) 
5  4  5  bo 

100. 4 2094b.- 125. 469 5b ,-104 5. 7642b.  =  0 


f 8  Cai *bi > 


+ 


b.,+  11  .  301 10 5b .+2  3. 009 17 6b. +110. 95452b, 

3  4  S  6 

63. 78098 (b.,-a-)-0.  34346a. -22 .875617a-21 1.89899a,  = 
If  4  5  o 


where  i  *  0,1,... ,7. 


Now,  substituting  the  values  of  and  in  terms  of  x^. 


i  =  l, 2,..., 7  from  (32. A)  yields  a  set  of  nonlinear  equations  shown  be¬ 
low.  It  is  noticed  that  as  a^  =  the  equation  =  ao~^0  =  ® 

gives  no  information.  The  rest  of  the  equations  are 

f  (x , ....x?)  =  x2x^+0.6318A22396[x^(x^-x^)-3.98319992x2x^]  =  0 

(33.1) 

f2(x  , . . . ,x7)  =  x^ (3 . 22822291x2+8. 522553136x^-6 .939879587x^ 

+  x,-l)+x0(l. 582676549x,-13. 17807554x, ) 

.5  2  o  -> 

+  x,  (x.-x., ) -3. 983199922 (x0+x. Xr)  =  0  (33.2) 

6  A  1  2  1b 

f3(x  . x?)  =  x2(-12.71361621x6-x5)+x7(13.58355291x1 

+  1. 820964317x„-26. 29716913x. )+10.798AA539(x1x,+x.) 

2  A  16  7 

-  13.31248704(x4x6+x3x7)+6.327224282(x]+x5) 

+  1.588477708x6(l-x3)+20.035271A3(x2+x1x5)  =  0 

(33.3) 

f . (x, , . . . , X- )  =  x  (x„+668 . 46 700 7 lx  -281 . 4809 Ax  )+x, ( 386 . 986067 3x9 

4  1772  4  16  2 

+  362. 767005-456. 258273x3)-647.2403649(x4x6+x3x7) 

-  57 . 53603068x2x5-548 .5188427  (x^x^) 

+  2 35. 861385 (x1x6+x7)+235. 2945185+590. 5096275 (x3 
+  x  )  -  0  (33.4) 

fr (x, ,. . . fxJ  =  2357. 408023(x.x,+x7)+x, (1598. 839931x,+17096. 15228 

5  1/  1  0  /  o  2. 

-  32881. 95043x3)+x7(4.16745091x2+1599.839931x1-x4) 

-  472 . 67 35 322 (x. x,+x_x7)+24996 . 98242-9 39.0287936 (x. 

U  b  j  /  i 

+  x5)-2765.323026(x2+x]x5)-52.07771943x2x5  =  0 


(33.5) 
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f  (x  , . . . ,x7)  =  x, (-99. 4209415x_+11132.91981x„-57256. 87822) 

b  1  /  b  l  1 

+  X- (x. -100 . 4209415x, )+411 . 4274907 (x . x.+x.x.) 

/  4  I  4  6  1/ 

+  67006. 93001  (x^x^+1234. 5674  33(x2+xx) 

+  42.69011171X.X,  +  2 3203. 53455-39112 . 69694 (x  x, 

/  j  lb 

+  x?)  =  0  (33.6) 

f7(x  ,...,x7)  =  496.854897(xox,+x1x,)+198512.9704(x.x,+x,) 

/  1  /  Z  0  1/  1  o  / 

+  2228495. 69 5xc-l 70. 649 7831 (x. x  +x.x  ) 

6  4  6  3  7 

-  77618. 59617x,x,-3442861. 087-395845. 4335 (x,+xc) 

lb  .15 

-  8390.812346(x2+x1x5)-62.08251489x2x5  =  0  (33.7) 

Equation  (33)  is  a  set  of  high  order  nonlinear  simultaneous  equations 
which  is  very  difficult  to  solve.  However,  with  proper  initial  guesses 
the  Newton-Raphson  [15]  method  can  be  applied  to  solve  it.  Therefore, 
the  problem  lies  in  finding  an  appropriate  set  of  initial  guesses.  In 
this  case,  the  following  method  is  suggested  for  estimating  the  initial 
guesses. 

As  mentioned  earlier,  the  block  diagram  of  the  structure  of  the 
desired  fixed  configuration  control  svstem  is  shown  in  Fig.  5-1.  With¬ 
out  affecting  the  overall  transfer  function  of  the  system,  this  struc¬ 
ture  can  be  modified  into  a  form  as  shown  in  Fig.  5-2.  The  overall 
transfer  function  of  this  structure  is 

T  (s)  =  T  (s)  p-—  (34-1} 

1  2  C2(s) 

where  * 

G_  (s)G  (s)G_(s) 

T  (s)  -  -± - -2- - 

1+G1(s)G2(s)G0(s) 


W 
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Where  Gq(s)  ,  the  approximate  transfer  function  of  Gq(s),  is 
given  in  (31.5). 

The  purpose  is  to  determine  G^(s)  and  G^(s)  such  that  the  re¬ 
sponse  of  T^(s)  becomes  close  to  that  of  the  standard  model  T^Cs)  in 
(9).  Replacing  the  series  compensator  G^(s)G^(s)  in  Fig.  (5-2)  by  the 

designed  stabilization  filter  F  „(s)  in  (30)  the  resulting  transfer 

sz 

function  T^(s)  in  (34.1)  is  equated  to  the  standard  model  Tr(s)  in  (9) 
as  follows. 


+175816571. ls2+425620,1128s3 
+ 205208030. 9s2+215915050. 5s3 

+ 154492. 684s4+29891. 09151s5 

- r~i] 

+  117. 470784s  +s 


(34.2) 
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Substituting  Tr(s)  in  (9)  into  (34.2)  and  simplifying,  the  appropriate 
* 

transfer  function  G2(s)  of  G2(s)  is  obtained. 

v  S.n36619203xl09+3.46495752xl010s+4.540060393xl010s2 
i s )  =  6  in  10  2 

Z  3.036619205x10  +3.008227329x10  8+4.613716606x10  s 

+7.840679235xl09s3+4.363076841xl09s4+1.763524302xl08s5 
+6. 124169121xl09s3+4.498497844xl09s4+8.512494768xl07s5 

+4.2562Q1128x105s6 _ : _ 

+9.985459768xl05s6+9.698650697xl03s7+49.1568119s8+0.243466s9 

(34.3) 

* 

A  set  of  dominant  quotients  h^  of  G^s),  given  below,  are  de¬ 
termined  by  expanding  (34.3)  into  a  continued  fraction  of  the  second 
Cauer  form 


h  =  -1.102755917 

h.}  =  -0.1287948973 

h,  =  5.593229805 

u 

h5  =  0.1338916858 


(34.4) 


Substituting  the  first  five  quotients  hj,h2,...,h5  into  (13.5) 

*  * 

gives  a  second  order  approximate  model  G2(s)  of  the  approximate  parallel 
filter  G2(s)  in  (34.3)  as 
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,  0.994929057s2+0.739A973923s2+0. 1058245527  _ 

G  (s)  =  — r - (34.5) 

s  +0.643533679s+0. 1058245527 

**  * 

(s),  the  approximate  model  of  G2(s),  is  also  an  approximate  model  of 

G,(s)  in  (32.2). 

* 

The  approximate  model  G^(s)  of  the  series  compensator  G^(s)  in 
Fig.  5-1  can  be  obtained  as  follows 


Fs2(s)  =  856 ■ 62 8 596s 3+2 1834 . 4682 ls2+l 3787 ,1076s 
G**(s)  0.994929057s4+4. 040722745s3 

+2252. 284999 _ 

+1 3237. 10512s2+9837.144919s+l407. 678125 


To  obtain  a  set  of  dominant  quotients  Equation  (35.1)  is  expanded  into 
a  continued  fraction  of  the  second  Cauer  form.  Some  of  the  quotients 
obtained  are 


h  =  0.625 

h2  =  1.845828612 

h3  =  0.0839039052 

h4  = 

h8  =  ...  (35.2) 


The  first  three  quotients  h  ,  h2 ,  h^  are  substituted  into  (13.3), 

A*  * 

which  gives  (s) ,  an  approximate  model  of  (s)  as  well  as  of  Gj(s) 
in  (32.1).  G^  (s)  obtained  is 


**  ,  ,  _  1. 41 0628426s+0. 2184671685 
1  "  s+0. 1365419803 
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(35.2) 


Comparing  (32.2)  with  (34.5)  and  (32.1)  with  (35.2)  we  have  a  set  of 
initial  guesses  to  solve  the  set  of  high  order  nonlinear  equations  in 
(33).  Thus,  the  set  of  initial  guesses  is 


and 


x*  =  0.643533679 

x*  =  0.1058245527 

x*  =  0.994929057 

x*  =  0.7394973923 
4 

x*  =  0.1 3» 54 19803 
x*  =  1.410628426 

O 

x*  =  0.2184671685 


(36) 


Using  these  initial  guesses  the  Newtcn-Raphson  method  [15]  is 
applied  to  solve  the  nonlinear  equations  in  (33).  It  is  found  that  the 
Newton-Raphson  method  converges  to  the  desired  solutions,  given  below, 
at  the  14th  iteration  with  the  error  tolerance  of  10  .  The  solutions 

are 


xL  =  0.503850 
x2  =  0.059928 
x3  =  1.051503 


x,  =  0.580016 
4 
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x5  =  4.831826 
x,  =  1.885577 

O 

x  7  =  6.744450 


Therefore,  the  desired  compensators  G^(s)  and  G^Cs)  are 


.  .  l,885577s+6. 744450  1 . 885577 (s+3. 57688) 

1(S;  s+4. 831826  "  s+4. 831826 


(37.1) 


and 


C,  (s)  = 


1.051503s  +0.580016s+0. 059928  1 . 051503(s+0. 13769) (s+0. 41391) 


s+0.503850s+0. 059928 


(s+0. 19244) (s+0. 311405) 
(37.2) 


The  unit  step  response  curves  of  the  existing  stabilized  system 
T^(s)  in  (1.5)  and  the  redesigned  system  using  the  compensators  G^(s) 
and  G^(s)  in  (37),  and  Gfl(s)  in  (1.3),  are  compared  in  Fig.  4.  The  re¬ 
sult  is  satisfactory. 

It  is  interesting  to  note  that  G^(s)  and  C^s)  in  Eq.  (37)  are 
positive  real  functions  with  positive  real  poles  and  zeros,  which  makes 
it  possible  to  realize  the  compensators  G^(s)  and  G^s)  using  passive 
elements,  whereas,  the  existing  stabilization  filter  a  non~ 


positive  real  function  and  it  is  realized  by  using  active  elements. 


CHAPTER  VI 


CONCLUSION 

The  existing  stabilized  pitch  control  system  has  been  redesign¬ 
ed  by  redesigning  the  existing  stabilization  filter.  Two  computer- 
oriented  methods,  a  dominant  data  matching  method  and  an  algebraic 
method,  have  been  presented  to  redesign  the  existing  stabilization 
filter.  Thus,  various  low-order  stabilization  filters  have  been  ob¬ 
tained.  As  a  result,  the  implementation  cost  of  the  missile  system  is 
reduced . 

The  proposed  dominant-data  matching  method  can  be  used  for 
various  purposes.  For  example,  when  the  specifications  or  the  design 
goals  of  a  control  system  are  given,  the  proposed  method  can  be  used 
to  obtain  a  standard  transfer  function,  which  significantly  simplifies 
the  design  process  in  the  frequency  domain.  When  a  high-order  transfer 
function  is  given,  various  low-order  models  can  be  obtained  with  the 
help  of  the  dominant-data  matching  method.  The  method  can  be  used  in 
the  problems  of  identification  as  well.  The  great  advantage  of  this 
method  is  that  the  transfer  functions  obtained  by  using  this  method  have 
the  exact  assigned  frequency-domain  specifications. 

The  algebraic  method  has  been  applied  to  achieve  the  advantages 
of  the  feedback  control  system  so  that  the  performances  of  the  redesign¬ 
ed  pitch  control  system  ran  be  greatly  improved. 

The  application  of  the  dominant-data  matching  method  always 


gives  rise  to  a  set  of  nonlinear  equations  which  can  be  solved  if  a  set 
of  proper  initial  guesses  is  known.  In  this  connection,  various  methods 


have  been  discussed  for  estimating  a  set  of  proper  initial  guesses. 

Finally,  it  is  important  to  mention  that  the  proposed  computer 
aided  design  methods  can  be  used  to  design  general  control  systems. 
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A  matrix  in  the  block  Schwarz  form  and  the  stability  of 
matrix  polynomials 

LEANG-SAN  SHIEHf  and  SHA1LENDRA  SACHETIf 


A  matrix  which  consists  of  block  elements  is  established  in  the  block  Schwarz  form 
via  a  linear  transformation.  The  transformation  matrix  constructed  by  Chen  and 
Chu  is  modified  anti  extended  for  transforming  the  block  companion  form  to  the 
block  Schwarz  form.  A  sufficient  condition  is  derived  for  determining  the  stability 
of  a  multivariable  system  whose  characteristics  arc  expressed  by  a  matrix  polynomial. 
The  matrix  polynomial  may  or  may  not  be  symmetric. 


1.  Introduction 

The  properties  and  applications  of  the  Schwarz,  matrix,  which  has  scalar 
elements,  has  been  investigated  by  various  authors  (Schwarz  1 956,  Parks  1 963, 
Wall  1948,  Anderson  et  al.  1976,  Barnett  and  Storey  1970),  and  the  transforma¬ 
tion  matrix,  which  relates  various  matrix  forms  and  the  Schwarz  form,  has 
also  been  established  by  numerous  authors  (Butchart  1965,  Chen  and  Chu 
1966,  1967,  Barnet  and  Storey  1967,  Loo  1968,  Power  1969,  Datta  1974, 
Sarma  et  al.  1968).  Chen  and  Chu  (1966,  1967)  constructed  a  transformation 
matrix  which  links  the  Schwarz  form  and  the  companion  form  by  using  the 
Routh  array  elements  (Routh  1877).  However,  all  existing  methods  (Schwarz 
1956,  Parks  1963,  Wall  1948,  Anderson  et  al.  1976,  Barnett  and  Store}-  1967, 
1970,  Butchart  1965,  Chen  and  Chu  1966,  1967,  Loo  1968,  Power  1969,  Datta 
1974,  Sarma  et  al.  1968)  deal  only  with  the  system  matrix  which  has  scalar 
elements  but  not  block  elements.  In  this  paper  a  matrix  which  consists  of 
block  elements  is  constructed  in  the  block  Schwarz  form  and  a  linear  transforma¬ 
tion  matrix  which  consists  of  block  elements  is  established  to  transform  the 
matrix  in  the  block  companion  form  (Shieh  1975)  (the  block  phase  variable 
form)  to  the  block  Schwarz  form.  A  sufficient  condition  is  then  derived  to 
determine  the  stability  of  a  multivariable  system  whose  characteristics  are 
described  by  a  matrix  polynomial  (Shieh  1975,  Shieh  et  al.  1976).  The  matrix 
polynomial  may  or  may  not  be  symmetric. 

2.  A  transformation  for  a  matrix  in  the  Schwarz  block  form 

Consider  a  set  of  linear,  time-invariant,  ordinary  differential  equations  in 
the  differential  matrix  polynomial  form 

"£  RiD'->A'(<)  =  (0],  Bn+l  =  /  (la) 

»=  1 

D*-'X( 0)  =  [«,_,],  i  =  1,  2,  rt  (16) 

where  X '(/)  is  the  m -dimensional  state  of  the  system.  The  B,-  are  m  x  m  real 
constant  matrices  and  the  differential  operator  D  is  D  =  dldt.  The  matrix  I 
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is  an  identity  matrix  and  [0]  is  a  null  matrix.  The  corresponding  state  equa¬ 
tion  of  eqn.  (1)  in  the  block  companion  form  is 

[*]  =  [£][*]  (2  a) 

WO)] -[at]  (2  b) 

where 


The  dimensions  of  the  matrix  [£],  each  block  element  in  the  [£],  and  the  state 
vector  [a-]  are  (n  x  m)  x  (n  x  m),  waxm,  and  (»xr»)xl,  respectively.  The  [£] 
is  the  matrix  in  the  block  companion  form  or  the  block  phase  variable  form 
(Shieh  1975). 

Equation  (2)  can  be  transformed  into  the  block  Schwarz  form  by  using  the 
following  linear  transformation  [K^]  : 

M-M*]  (3) 

and 

tfl-WilWlWiJ-’M 

-Wily]  (4) 

where 
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and 


[A]- 


0 

-A, 

0 

0 

0 


/  0  0 
0  7  0 

-A,  0  l 
0  0  0 
0  0  0 


0 

0 

0 

0 

'  An_  i 


0 

0 

0 

I 

-A„ 


The  dimension  of  each  block  element  in  the  matrix  [A }  and  the  matrix  [Arj]  is 
mxm.  The  [^1]  is  the  matrix  in  the  block  Schwarz  form.  The  linear  trans¬ 
formation  matrix  [A',],  which  relates  the  coordinates  [a:]  and  [?/]  in  eqns.  (2) 
and  (4),  is  constructed  by  following  the  approach  proposed  by  Chen  and  Chu 
(I960).  The  block  elements  C(j  having  dimension  mxm  in  eqn.  (3)  can  be 
obtained  from  the  following  matrix  Routh  algorithm  and  the  matrix  Routh 
array  (Shieh  and  Gaudiano  1974,  Shieh  et  al.  1976). 

Before  performing  the  matrix  Routh  array  we  define  l  =  (n/ 2)+l  if  n  is 
an  even  number  ;  otherwise,  Z  =  m  +  1/2,  and  the  double  subscripted  block 
elements  Cl>j  and  C2  j  become  : 

C\,j  —  &n+3-2 j>  j  —  1  >  2,  3,  . . . ,  l  | 

=  j =  1 ,  2,  3,  . . .,  I  >  (5  a) 


The  Bt  are  the  mxm  real  constant  matrices  shown  in  e  pi.  (1).  The  matrix 
Routh  array  and  the  matrix  Routh  algorithm  are 


h>=c21c31-' 

H3  =  C3]Cix  1 

L Cc 11 
/  °21 

yc3l 

ycu 

^12 

^22 

^32 

^42 

c 

C13 

c23 

C33 

^43 

- 

C24  . . . 

Ht-CuCn  1 

/°51 

°52 

c 

He  =  CgiC?!-1 

/  °61 
z<?71 

°82 

• 

• 

H  n  —  Cn,lCn+l, 

-1  /®n,  1 

^»+ 1 

1 

. 

where 

@i,j  —  —  Hi-zCi-ij+v  j=  R2,...,  t'*3,  4,  ... 

=  i(C<+i,i)_1>  i=I 

det  (Cf+ i  j)  #  0 


(5  b) 


The  matrices  H{  in  eqn.  (5  b)  are  called  the  matrix  quotients.  Performing  a 
new  linear  transformation 


(•> 
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on  eqn.  (4)  yields  an  alternative  matrix  [/’]  in  the  block  Schwarz  form  as 
follows  : 


[z]  =  [£,][  A  ][tf2]-*[Z] 

-[^IM 


where 


[K  2]  = 


Cn,x  0 

d  ^n- 1,1 


0 

0 


and 


[f]  = 


0  Hn_x-' 

-H-'  0 


tfn-r1 


-^n-r1 

0 

o 

o 

o 


o 

o 

o 

o 


o 

o 


o 

o 

0 


0  0 
0  0 


C21  0 


0  c. 


0  0 
0  0 
0  0 


■H-  4* 


H 3-‘ 
0 

-h3-' 

0 


0 

H 


(7  a) 


(7  6) 


0 


Hr1 


(7  c) 


The  [if2]  is  a  block  diagonal  matrix  having  the  diagonal  block  elements  ob¬ 
tained  from  the  block  elements  C{1,  i=  1,  2 . which  are  in  the  first  column 

of  the  matrix  Routh  array  in  eqn.  (5  6),  while  the  matrix  [P]  is  the  required 
matrix  in  the  block  Schwarz  form  which  can  be  constructed  by  using  the  matrix 

quotients  Ht,i=  1,2 .  obtained  from  the  same  matrix  Routh  array.  A 

similar  matrix  (Schwarz  and  Friedland  1965),  which  was  formulated  in  the 
Schwarz  form  but  not  in  the  block  Schwarz  form,  was  used  to  prove  the 
stability  of  a  linear  system  by  Parks  (1963). 

The  linear  transformation  matrix  [A],  which  links  the  coordinates  [x]  in  the 
block  companion  form  and  the  coordinates  [z]  in  the  block  Schwarz  form,  is 


(8  a) 


=  [K][x] 
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where 


IH  1. 1 

0 

0 

0 

0 

0 

0 

H  H  /  n  1,2 

0 

0 

0 

0 

0 

0 

0  •  • 

0 

0 

0 

0 

0 

"  n  ■/  n  3,3 

0 

H  rJ  61 

0 

0 

0 

0 

0  •  • 

52 

0 

U  4^51 

0 

0 

0 

//  /»  .  . 

11  n  -6*  n  -5.4 

0 

42 

0 

H/  41 

0 

0 

o  •  • 

33 

0 

II./  32 

0  ! 

H  31 

0 

0 

0 

22  i 

0 

"A,  . 

The  matrix  [A']  is  a  block  triangular  matrix.  All  the  block  elements  in  eqn. 
(8  b)  can  be  directly  obtained  from  the  matrix  Routh  array  in  eqn.  (5  b).  For 
example,  the  block  elements  in  the  main  diagonal,  which  are  shown  by  the  first 
dotted  diagonal  line,  are  obtained  by  the  respective  products  of  the  matrix 
quotients  and  block  elements  C{1  in  the  first  column  of  the  matrix  Routh 
array.  The  block  elements  of  the  first  lower  diagonal  in  the  [A]  are  null 
matrices,  and  the  block  elements  of  the  second  lower  diagonal  in  the  [A],  which 
are  shown  by  the  second  dotted  diagonal  line,  are  determined  by  the  respective 
products  of  the  matrix  quotients  Hi  and  the  block  elements  Ci  t  in  the  second 
column  of  the  same  matrix  Routh  array,  etc.  The  sizes  of  the  matrices  [A] 
and  [A]  are  determined  by  the  degree  of  the  matrix  polynomial  and  the  order 
of  the  matrix  coefficients  in  eqn.  (1).  For  instance,  when  the  degree  of  a 
matrix  polynomial  is  »  =  4  and  the  dimension  of  each  matrix  coefficient  is  m, 
the  corresponding  4mx4m  square  matrices  [A]  and  [A]  are  taken  from  the 
right-hand  side  lower  corner  of  the  matrices  [A]  and  [A]  in  eqns.  (7  c)  and  (8  b). 


3.  A  sufficient  condition  for  the  stability  of  a  matrix  polynomial 

In  a  single  variable  system  the  Routh  criterion  (Routh  1877)  is  applied  to 
the  characteristic  polynomial  of  a  linear  system  for  determining  the  absolute 
stability.  In  other  words,  the  scalar  polynomial  in  the  form  of  eqn.  (1)  is 
arranged  in  the  Routh  array  in  eqn.  (5  b),  then  the  Routh  criterion  is  applied. 
If  the  scalars  C( ,  in  the  first  column  of  the  Routh  array  have  no  sign  changes  or 
all  elements  Ci  v  i=  1,  2,  . ..,  are  positive  real,  then  the  system  is  asymptotically 
stable.  Since  the  Routh  algorithm  and  the  Routh  array  have  been  successfully 
extended  to  the  matrix  Routh  algorithm  and  the  matrix  Routh  array  (Shieh 
and  Gaudiano  1974,  Shieh  et  al.  1976),  also  a  positive  definite  matrix  (Bellman 
1970)  is  commonly  considered  as  a  natural  generalization  of  a  positive  number. 
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it  is  interesting  to  ask  whether  or  not  a  multivariable  system  whose  charac¬ 
teristic  matrix  polynomial  shown  in  eqn.  (1)  is  asymptotically  stable  if  the 
block  elements  C(>1,  t=  1,  2,  in  the  first  column  of  the  matrix  Routh  array 
in  eqn.  (5  b)  are  positive  definite  matrices.  In  other  words,  can  we  directly 
extend  the  Routh  criterion  (Routh  1877)  to  the  matrix  Routh  criterion  ? 
This  paper  will  partially  answer  this  question. 

Because  the  stability  of  a  system  is  invariant  under  the  transpose  operation 
of  the  system  matrix,  we  consider  the  following  system  : 

[q]  =  [F1][q] 

«[G]fo]  (»> 

The  matrix  [F]  in  eqn.  (9)  is  defined  in  eqn.  (7)  and  the  transpose  of  the  matrix 
[F]  is  defined  as  [G].  If  the  matrix  quotients  Ht  in  eqn.  (5  b)  are  positive- 
definite  symmetric  and  real  matrices,  then  we  can  consider  the  following 
quadratic  equation  (Kalman  and  Bertram  I960,  Barnett  1971)  : 

f=[<zTl[P]fo]  (10 «) 

where 

~Hn  0  -00 

0  Hn •  0  0 

[P}= . 

0  0  •  Ht  0 

0  0  •  0  Hl 

The  derivative  function  v  is 

v  =  [qt][PG+GrP][q] 

=  -[?TIQ][?]  (10  b) 

where 

'0  -I  0  •  0  0  I  [0  0  0  •  0  0 

1  0  -/  •  0  0  000-00 

[P][G]=  0  1  0-00,  [Q]=  0  0  0-0  0 

00  0-0  -I  000-00 

0  0  0  •  I  -l\  |_0  0  0  •  0  2/ 

The  v  function  in  eqn.  (10  a)  is  in  a  positive -definite  quadratic  form  and  the  v 
function  in  eqn.  (10  b)  is  in  a  negative -semidefinite  form.  Therefore  the  system 
in  eqn.  (9)  or  in  eqn.  (1)  is  asymptotically  stable.  From  the  above  derivation 
we  obtain  the  sufficient  condition  that  a  multivariable  system,  whose  charac¬ 
teristic  matrix  polynomial  has  the  form  shown  in  eqn.  (1),  is  asymptotically 
stable  if  the  matrix  quotients  H(  in  eqn.  (5  b)  are  real  symmetric  positive- 
definite  matrices.  From  eqns.  (2  a)  and  (7  c)  it  can  be  observed  that  the 
Bn(  =  H1~l  =  Cn)  must  be  symmetric  and  positive-definite  for  the  sufficient 
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condition  to  apply.  It  is  also  noted  that  this  sufficient  condition  deals  only 
with  //;  and  not  C\  j.  This  implies  that,  if  a  system  is  asymptotically  stable, 
the  block  elements  Ciml,  »=1,  2,  ...,  in  the  first  column  of  the  matrix  Routh 
array  and  the  Bt  in  eqn.  (1)  are  not  necessarily  symmetric  and  positive- 
definite  matrices.  In  other  words,  a  positive-definite  matrix  is  not  necessarily 
a  natural  generalization  of  a  positive  number,  and  the  necessary  and  sufficient 
condition  of  the  Routh  criterion  (Routh  1877)  cannot  be  completely  extended 
to  the  matrix  Routh  criterion  for  general  eases. 

On  the  other  hand,  the  necessary  conditions  for  asymptotic  stability  due 
to  the  Routh  criterion  (Routh  1877)  can  be  partially  extended  to  the  case  of 
matrix  polynomials.  The  necessary  conditions  are  described  as  follows  : 

(i)  The  determinant  of  B1  in  eqn.  (1)  is  non-zero. 

(ii)  The  determinants  of  Bn+1  and  Bx  in  eqn.  (1)  have  the  same  sign  if  the 
determinant  of  Bn+1  (  =  CU)  is  non-zero. 

These  conditions  can  be  easily  verified  from  the  basic  laws  of  algebra.  Thus, 
in  this  paper,  we  have  partially  extended  the  Routh  criterion  (Routh  1877)  to 
the  matrix  Routh  criterion  for  determining  the  asymptotic  stability  of  a  class 
of  matrix  polynomials. 

Sometimes  in  applying  the  approach  proposed  in  this  paper  difficulties  may 
be  encountered  in  calculating  the  matrix  quotients  Ht  in  eqn.  (5  b).  This 
implies  that  if  any  Ci  t  in  eqn.  (5  b)  is  singular,  then  the  Hi  cannot  be  obtained 
to  determine  the  stability  of  a  matrix  polynomial.  This  limitation  can  be 
remedied  by  multiplying  the  original  matrix  polynomial,  defined  as  T(s),  by  a 
polynomial  matrix  defined  as  E(s),  where  the  roots  of  the  determinant  E(s)  have 
negative  real  parts.  Then  we  apply  the  matrix  Routh  procedure  to  the 
modified  matrix  polynomial  T{s)E(s)  or  E(s)T(s).  It  is  noted  that  the  stability 
of  the  original  system  is  reserved  because  the  stability  of  a  system  is  invariant 
under  this  modification.  An  alternative  way  is  to  perform  transformation 
s— >l/s  to  the  T(s)  and  then  applying  the  matrix  Routh  procedure  to  the 
modified  matrix  polynomial  defined  as  Tj(s)(=  T(s)|s_,/S).  In  other  words, 
the  C'u  and  C2j  in  eqn.  (5  a)  are  rewritten  as  follows  : 

('ij  =  B2j_i  for  j  =  I,  2,  3,  ... 

,}  —  for  j  =  1 ,  2,  3,  ... 

Again,  the  stability  of  the  original  system  is  invariant  to  this  modification  and 
the  numerically  ill-conditioned  problem  (i.e.  C,-  j  is  singular)  can  be  solved. 
Examples  are  established  in  this  paper  to  show  the  procedure. 

4.  Illustrative  examples 

4.1.  Example  1 

Consider  the  following  differential  matrix  polynomial : 

"l" *  BtDi-'X(t)  =  [Q]  (11) 

•  ■  i 
or 

(«  (9) 

B3X(t)  +  BtX{t)  +  B3X(t)  +  B2X(t)  +  B2X(t) 

(4)  (S> 

=  CnX(t)  +  C21X(t)  +  C12X(t)  +  C22X(t)  +  C13X(t)  =  [0] 
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where 


A  matrix  in  the  block  Schwarz  form  and  the  linear  transformation  matrix  which 
transforms  the  block  companion  form  to  the  block  Schwarz  form  are  required 
to  be  constructed,  and  the  stability  of  the  system  is  of  interest. 

Arranging  the  matrices  R,  in  eqn.  (11)  in  the  matrix  Routh  array  in  eqn. 
(5  b)  results  in  the  following  : 


The  matrix  quotients  H{,i=  1,2 . 4,  in  eqn.  (12)  are  positive-definite 

symmetric  real  matrices  ;  therefore  the  system  is  asymptotically  stable.  It  is 
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noted  that  the  block  elements  C{1,  i  =  1 . 5  in  the  first  column  of  the  matrix 

Routh  array  in  eqn.  (12)  are  not  all  positive-definite  symmetric  real  matrices. 
The  state  equation  in  the  block  companion  form  is 

[i]  =  [£][*]  (13) 

where 


and  the  state  equation  in  the  block  Schwarz  form  is 

[i]  =  [F][z]  (14  a) 

where 


It  is  interesting  to  note  that  the  characteristic  equation 


\el  -5|  =  |s/-Jf|=0 
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and  the  roots  which  have  negative  real  parts  are,  respectively  : 

s8  4-  3s7  +  28-95s«  +  79-35s6  +  206s4  4-  458-875 s3 

4- 22  l-05s2  4-48-4*4- 4  =  0  (146) 

and 


-0-0239155  ±j4-27199 
-  0-0784809  ±j2-95637 

-0-189163  ± jO-165319  (14  c) 

-0-177194 

-2-23969 

The  linear  transformation  matrix  between  the  [ar]  and  [z]  coordinates  is 

[z]  =  [K][x]  (15) 

where 


4.2.  Example  2 

Consider  the  following  transfer-function  matrix  [71(«)]  which  is  expressed 
by  the  product  of  two  relatively  prime  polynomial  matrices  [2V(s)]  and  [Z)(s)]_1 
or 

[T(s)]  =  [#(*)][!)(*)]-»  (16) 

The  characteristic  matrix  polynomial  [Z>(s)]  is 

[  D(s )  ]  =  Bhs*  +  Bis3  +  Bts2  +  Bj$  +  Bx 

—  4“  C21S®  4"  4-  4"  Cjg 
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where 


It  is  interesting  to  determine  the  stability  of  this  system.  Following  eqn. 
(5  6)  yields  the  matrix  Routh  array  as  follows  : 


The  matrix  quotients  H{,  i—  1,  4,  in  the  matrix  Routh  array  are  positive- 

definite  symmetric  real  matrices.  Therefore  the  system  is  stable.  It  is 
observed  that  the  block  element  £31  in  the  first  column  of  the  matrix  Routh 
array  is  not  symmetric. 


256 


L.  S.  Skieh  and  S.  Sacheti 


4.3.  Example  3 

Consider  the  stability  and  the  structure  of  the  matrix  Routh  array  of  the 
following  matrix  polynomial  T(s)  are  of  interest  : 

T(s)  =  B^s3  4  B^s3  4  B^s  +  B1 

=  Cn' s3  4  C21'  s2  4  C12' s  4  C22'  =  0  (18a) 

where 


The  determinant  of  Bt(  =  Cn')  =  -l  and  that  of  i?,(=C22')=  1.  From  the 
derived  necessary  conditions  for  asymptotic  stability  we  conclude  that  the 
system  is  unstable  because  the  determinants  of  Bt  and  Bl  have  different  sign. 
Further  study  of  the  stability  is  not  necessary.  It  might  be  interesting  to  see 
the  corresponding  characteristic  equation  of  this  system  which  can  be  written 
as  follows  : 


det  T(s)=  -s«  — 2s5 4 3s2 4 2s 4  1=0  (18  6) 

Because  the  first  and  the  last  coefficients,  which  are  the  determinants  of  Bt 
and  Bl  respectively,  have  different  sign,  therefore  the  system  is  unstable. 
In  order  to  study  the  structure  of  the  matrix  Routh  array  of  this  unstable 
system  and  the  numerically  ill-conditioned  problem  (i.e.  Ci  x  is  singular)  we 
apply  the  matrix  Routh  algorithm  in  eqn.  (5)  and  use  the  C,-  /  in  eqn.  (18  a). 
The  matrix  Routh  array  cannot  be  obtained  because  C21'  is  singular.  This  is  a 
numerically  ill-conditioned  case.  Since  the  stability  is  invariant  between  the 
original  system  T(s)  and  the  modified  system  T1{s)(  =  T(s)\,^1is),  we  can 
construct  the  matrix  Routh  array  for  the  Tj(s).  The  T,(s)  can  be  written  as 
follows  : 


where 


Ti($)  —  T(s) |s_ifo—  C*22  334C12  *24C21  S4C,, 
=  Cns3  4  C>2  4CjjS4C22  =  0 


cn  = 


f 

C,2=( 

(' 

') 

\° 

1/ 

(18c) 
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The  corresponding  matrix  Routh  array  is 


X  0 

h, =cuc2r1= 

\o  1 


/I  0 

tf*  =  C2I<Vl  = 

\o  1 


/O  1 

H3  =  C3lCir'  =  [ 

\l  0 


Although  the  C12  is  singular,  we  can  determine  all  the  s.  It  is  observed 
that  the  //,  and  H2  are  symmetric  and  positive  definite  matrices,  while  the 
H3  is  a  symmetric  and  non-positive  definite  matrix. 

An  alternative  method  can  be  described  as  follows.  Let  us  construct  a  new 
matrix  polynomial  T2(s)  by  multiplying  a  matrix  polynomial  E(s)  =  (s  +  1)/  to 
the  T(s)  and  then  defining  the  matrix  coefficients  as  Chi'  and  C2i'  : 

T2(s)  =  (a  +  1 ) T(s)  =  Cu' s*  +  C21'  s3  +  C12  s 3  +  C22'  s  +  C13'  =  0  (18  e) 

where 


If  we  wish  to  maintain  the  consistency  of  Cn  =  /,  we  may  interchange  the  rows 
in  the  T2(s)  and  define  new  matrix  coefficients  as  Chi  and  C2  i : 


T2'(s)  =  Cns*  +  Ctls 3  +  Clts*  +  Csrs  +  Cl3  -  0 


(18/) 
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where 


e\,= 


The  corresponding  matrix  Routh  array  is 
C 

t  2  -l\j 

=  \ \ 


Hx-i 


//4=i 


\  /I 

2\ 

(° 

)  c’“‘(2 

1 

l) 

I  C„-| 

[l 

J 

Cu- 


/0  1' 


J  0 


6  -1; 


(18g) 


No  singular  matrix  appears  in  the  matrix  Routh  array  and  all  the  i/,’8  can  be 
obtained.  It  is  observed  that  only  the  H3  is  a  symmetric  but  non-positive 
definite  matrix. 

From  the  above  illustrations  we  conclude  that  if  any  ill-conditioned  problem 
occurs  in  the  calculation,  then  the  above  methods  can  be  applied  to  solve  the 
problem. 


5.  Conclusion 

The  transformation  matrix  established  by  Chen  and  Chu  (1966)  for  trans¬ 
forming  the  companion  form  to  the  Schwarz  form  has  been  modified  and 
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extended  to  transform  the  companion  block  form  to  the  block  Schwarz  form. 
The  new  matrix  in  the  block  Schwarz  form  has  been  constructed  by  using  the 
matrix  quotients  obtained  from  the  matrix  Routh  array  which  is  constructed 
from  the  characteristic  matrix  polynomial.  When  the  matrix  quotients  in  the 
matrix  Routh  array  arc  positive-definite  symmetric  real  matrices,  the  suffi¬ 
cient  condition  derived  in  this  paper  shows  that  the  multivariable  system  is 
asymptotically  stable.  Also,  a  set  of  necessary  conditions  has  been  derived 
for  the  asymptotic  stability.  Thus,  we  have  partially  extended  the  Routh 
criterion  (Routh  1877)  to  the  matrix  Routh  criterion  for  a  class  of  matrix 
polynomials.  The  direct  extension  of  the  necessary  and  sufficient  condition 
of  the  Routh  criterion  (Routh  1877)  to  a  genera)  matrix  polynomials  need 
further  studies. 
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Abstract— A  simple  method  is  proposed  that  will  fit  the  coefficients  of  a  transfer  function  from  the  real  and 
imaginary  parts  of  experimental  frequency-response  data.  An  approximate  logarithmic  amplitude-frequency 
plot  is  used  to  formulate  an  irrational  transfer  function  which  then  estimates  the  interpolation  data  and  the 
degree  of  the  final  transfer  function.  The  present  method  is  applicable  to  either  minimum  or  non-minimum 
phase  system  identification. 


I.  INTRODUCTION 

The  problem  of  finding  unknown  coefficients  of  a  transfer  function  as  a  ratio  of  two  frequency- 
dependent  polynomials  has  been  investigated  by  Levy[l],  Kardashov  and  Karniushin[2],  and 
Sanathanan  and  Koerner[3],  In  general,  they  would  evaluate  the  polynomial  coefficients  by 
minimizing  the  weighted  sum  of  squares  of  the  errors  in  magnitude  at  arbitrary  experimental 
points.  Ausman[4]  proposed  a  graphical  method  to  rapidly  estimate  the  coefficients  of  a  transfer 
function;  however,  that  procedure  is  only  applicable  for  a  minimum  phase  system. 

In  this  paper  a  simple  method  is  presented  to  approximate  the  coefficients  of  a  transfer 
function  for  minimum  and  non-minimum  phase  systems.  The  generalized  Bode  plot  is  used  to 
formulate  an  irrational  transfer  function  from  which  we  obtain  interpolation  frequency- 
response  data  that  will  allow  us  to  estimate  the  polynomial  coefficients  without  minimizing  the 
weighted  sum  of  squares  of  the  errors  in  magnitude  at  arbitrary  points. 

2.  THE  DERIVATION 

Consider  the  transfer  function 


ru  ,  _  Po  +  p,s+p2s2+---  +  pmsm 
I  +  q,s  +  q2s2  + - H  q„s” 


(1) 


where  p,  and  q,  are  unknown  coefficients  to  be  determined.  Substituting  s  =  j<ok  we  have 

„  .  _  (Po~  P2(ilk  +  P*Vl>*  ~  Pt,Wk  +  ~  ’  ')  +  j(PtO>k  ~  p-jtlik'  +  P;«Jfc5  -  Pr<l>k  +  ■  •  •) 

*  ( 1  -  q2a>k  +  -  qt<ok  +  •  •  ■)  +  j(qitok  -  q2(ok  +  qs<ok  -  q?^*  +  •  ■  •) 

=  R(<ok)  +  jI{tok) 

-  Rk  +  ih 


where  Rk  and  /*  are  the  real  and  imaginary  parts  of  the  transfer  function  at  the  experimental 
frequencies  mk.  After  we  multiply  both  sides  of  eqn  (2)  by  the  common  denominator,  we 
separate  the  real  and  imaginary  parts  and  then  equate  the  respective  real  and  imaginary  parts. 
We  now  have 

Po _ p2t»k  +  P*o>k  -  Pita*  +  •  •  ■  +  q ili&th  +  q-iRifiik  -  qj/kW -  qtRtfo, *  +•■  =  Rk  (3) 
and 

Pi“*  -  Pi***5  +  PfO>k  +  •  •  ■  -  q^Rucoh  +  q2Ihfok  +  qj/?*W  -  qJk<a*  +•■•  =  /».  (4) 
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Combining  eqns  (3)  and  (4)  results  in 


Po  +  Pi«*  -  Pitot1  -  Pyotk  +  p<<u  *  +  PiO>k - -  -  4i (Rk  -  h)tuk  +  q^Rk  +  Ik)to2  +  Qi(Rk  ~  /*)<o*5 

-  q*(Rk  +  /*)a>*4  =  Rk  +  Ik  (5) 


The  complete  form  of  eqn  (3)  is 


I  2  3  4  5 

1  OI|  ~  <Ol  ~  0)\  U>l  OJ\  . 

I  ah  0)2“  ~  ah^aii^'  • 

•  2  3  4  3 

lo>„-  U)m  -  0)m  U)m  (l)m  . 


~  T\o)i  S(W|"  Tinji3  —  svoj/  —  T\o>\  S|0>|6  .  - . 

—  TlO>2  S20>22  Tj< i>23  —  S2IO2*  ~  T 2(1)2  S20>2>  ■  ■  ■ 

Tm(OmSm(Om~7'm(Om  ~  Sm(Om  Tm(l)m  Sm(Om  •  •  • 


L  I  <D, 


■  0>k  -  0>k  <0*  (Ox 


~  TxOJx  SxO>x2  TxO)*  -  SxOt*  -  Tx(Ox  S^x 


where 


Po' 

'  R,  +  /, 

P\ 

R2  +  I2 

Pm 

= 

Rm  +  An 

«•  - 

Rx+1, 

(6) 


Sk=Rk  +  h\  k  =  1,2,... 

Tk  =  Rk-h\  *  =  1,2,... 

x =m +«+l 

By  substituting  the  selected  x  sets  of  frequency  response  data  into  eqn  (6),  we  can  solve  for  the 
required  unknown  coefficients  pt  and  q,. 


3.  ESTIMATION  OF  THE  CORNER  FREQUENCY  AND  THE  ORDER 
Bode  [3]  uses  piecewise  linear  segments  with  integer  slopes  to  approximate  the  logarithmic 
amplitude-frequency  characteristic  of  a  function.  Ausman[4]  applies  this  characteristic  to 
evaluate  the  coefficients  of  a  transfer  function.  Polonnikov[6,7]  generalizes  Bode's  approach  to 
estimate  the  phase-frequency  characteristic.  We  shall  now  obtain  a  logarithmic  amplitude- 
frequency  characteristic  by  piecewise  linear  segments  with  accurate  integer  or  fractional  slopes. 
The  approximate  transfer  function  is 


F(s)  = 


»o+a>fcr-('+r 
"  KrKr-Kr 


(7) 


where  at  and  bt  are  corner  frequencies,  and  where  nu  and  n,  may  be  integer  or  fractional 
values.  In  general,  eqn  (7)  is  an  irrational  transfer  function.  Compared  to  an  approximation 
made  by  other  methods  [4, 3],  this  present  analysis  is  much  better  because  the  slopes  may  be 
precise  fractional  values.  However,  the  worst  errors  caused  by  piecewise  segment  ap¬ 
proximation  occur  at  the  corner  frequencies  ax  and  bx\  therefore,  these  corner  frequencies 
provide  the  most  important  information  of  the  frequency-response  curve.  If  the  interpolation 
data  in  eqn  (6)  include  these  important  comer  frequencies,  a  good  transfer-function  fitting  is 
expected.  In  this  paper  these  corner  frequency-response  data  are  chosen  as  main  interpolation 
points  for  determining  the  unknown  coefficients  in  eqn  (6).  The  difference  of  the  order  of  two 
polynomials  in  eqn  (1)  can  be  estimated  from  eqn  (7).  In  other  words 


n-m 


(8) 


Based  on  eqn  (8),  the  numbers  of  the  unknown  coefficients  and  the  interpolation  points  may  be 
estimated. 
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4.  ILLUSTRATIVE  EXAMPLES 
Example  1.  Consider  Levy's  non-minimum  phase  example.  The  frequency-response  data 
generated  from  the  transfer  function  in  eqn  (9)  is  shown  in  Table  1  and  the  log-amplitude  plot 
versus  log-frequency  is  shown  in  Fig.  I. 


T(s)  = 


l-s 

l+0.ts  +  0.01sJ' 


(9) 


The  irrational  transfer  function  approximated  from  the  generalized  Bode  plot  is 


T(s)  = 


(■♦& 

054  ✓  „\0.57  /  V 

1  H)  (1+®) 

0.15 

1 

( 

'1+^ 

,  10J 

TB - 

where  the  corner  frequencies  are 


(10) 


tt|  =  «3=0.5,  0)3  =  <09  =  10 

<02=106  =  2  ,  <04=  toil  =  40 


Table  I. 


k 

<a'k 

1 

Given  data 

Identified  results 

|T(j»i)| 

/T( j«i) 

ft, 

|B| 

ft 

1 

0.1 

mm 

-6.45 

-0.1130 

1.0013 

-6.26 

0.9953 

-0.1092 

2 

msi 

1.0239 

-12.41 

mEJJI 

1.0160 

-12.41 

0.9923 

-0.2183 

3 

HS 

1.1194 

-29.43 

0.975 

musm 

1.1142 

-29.34 

0.9713 

-0.5459 

4 

0.7 

1.2393 

-39.01 

1.2171 

-38.91 

0.9472 

-0.7644 

5 

1.0 

1.4399 

-51.06 

-1.1200 

1.4125 

-50.66 

0.8955 

-1.0924 

6 

2.2772 

-75.04 

iEuki.ufl 

2.2631 

-75.15 

0.5798 

-2.1875 

7 

4.0 

4.4375 

-0.925 

4.3954 

-101.50 

-0.877 

-4.3071 

8 

7.0 

8.1751 

-135.9 

-5.870 

-5.6900 

8.0864 

-136.08 

-5.826 

-5.608 

9 

10.0 

-174.0 

-10.00 

-1.1050 

9.9115 

-174.67 

-9.869 

-0.9206 

10 

20.0 

5.5541 

-233.4 

-3.310 

4.460 

5.4612 

-233.4 

-3.245 

4.3926 

11 

40.0 

2.5451 

-253.5 

-0.7240 

2.4400 

2.5363 

-253.5 

-0.714 

2.4338 

12 

70.0 

1.4479 

-261.0 

-0.2270 

1.4300 

1.4205 

-261.0 

-0.225 

1.4026 

13 

100 

0.9994 

-263.5 

-0.1130 

0.9930 

0.9892 

.. 

-263.6 

-0.109 

0.9832 

20Io9,„|F(/»)|.  d 8 

20  loOiol<30">l-  dfl 


Tht  Identified  function:  fl(j>  -  rrti^?n  » 

t  ♦  010053s  +  0.010072*’ 

Fig.  I.  Bode  plot  shows  magnitude/frequency  response  and  piecewise  segment  approximations  of  F(j)« 

(l-»)/0+0.U+0.tsI). 
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The  order  of  eqn  (I)  may  be  estimated  from  eqn  (10),  or 

m  =0.54  +  0.57  +  0.15=1 
n  =2.25  =  2 
n  -  m  =  1. 


Four  frequency-response  data  (<■»■,  o>2,  m,  u*)  are  required  in  eqn  (6)  to  fit  the  four  unknown 
coefficients  p0,  Pu  4i  and  q2.  The  identified  transfer  function  is 


Ci  \  0.9%28-0,99l402i 

1  +  0.10053s  +0.0 10072s2' 


01) 


The  corresponding  frequency-response  data  of  eqn  (11)  and  that  of  eqn  (9)  are  compared  in 
Table  1.  The  results  are  very  satisfactory. 

Example  2.  A  set  of  frequency-response  data  generated  by  the  following  transfer  function 
is  shown  in  Table  2  and  the  log-amplitude  versus  log-frequency  plot  is  shown  in  Fig.  2. 


Table  2. 


*  «; 

Given  data 

Identified  results 

l  . 

/Tljmi) 

ft 

U 

K» 

it 

i 

0.1 

1.0002 

-0.2* 

1.0002 

-0.0048 

1.0002 

-0.28 

1.0002 

-0.0048 

2 

0.4 

1.0029 

-1.10 

1.0027 

-0.0193 

1.0029 

-1.10 

1.0028 

-0.0193 

3 

0.8 

1.0117 

-2.13 

1.0110 

-0.0375 

1.0124 

-2.18 

1.0116 

-0.0385 

4 

2.0 

1.1113 

-5.68 

1.1058 

-0.1 101 

1.1113 

-5.69 

1 . 1058 

-0.1 101 

3 

2.2 

1.1470 

-6.61 

1.1394 

-0.1321 

1.1470 

-6.62 

1.1394 

-0.1322 

6 

3.6 

1.4936 

-33.8 

1.2418 

-0.8299 

1.4935 

-33.8 

1.2416 

-0.8300 

7 

5.4 

0.8425 

-57.8 

0.4485 

-0.7132 

0.8424 

-57.8 

0.4484 

-0.7132 

8 

8.0 

0.6123 

-59.1 

0.3147 

-0.5253 

0.6123 

-59.1 

0.3147 

-0.5252 

9 

16 

0.3730 

-69.5 

0.1309 

-0.3493 

0.3730 

-69.5 

0.1309 

-0.3493 

10 

20 

0.3091 

-72.9 

0.0908 

-0.2955 

0.3091 

-72.9 

0.0908 

-0.2955 

II 

100 

0.0662 

-86.3 

0.0042 

-0.0661 

0.0662 

-86.3 

0.0042 

-0.0661 

12 

110 

0.0602 

-86.7 

0.0035 

-0.0601 

0.0602 

-86.7 

0.0035 

-0.0601 

20  log,0|G(/“)l.  «B 
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6.6378918s2  +  22.999878*  +  1 1 1.27974 
1  ”  s5  +  9.882741i J  +  28.37056*  + 111.27974- 


The  irrational  transfer  function  approximated  from  the  generalized  Bode  plot  is 


('♦d 

0.09 

1  1 

(■♦a 

^o.« 

(■♦d 

^  1.02 

(1+b) 

i"i 

H) 

"i 

<  s  ' 

,l+l06; 

TEW- 

1 

The  corner  frequencies  are 


«i  =  »ii  =  0.8,  ft)3  =  a>6  =  3.6, 

o>2  =  n»s =  2.2,  o>4  =  (os  =  5.4, 


Uf  —  <o\o—  16 
<0«  =  <0(2  ~  100. 


The  order  of  Eqn  (1)  is  estimated  as  follows: 


m  =0.09  +  0.63+1.02  =  2 
n=  2.42  +  0.28  +  0.04  =  3 
n  -m  =  1. 


(12) 


(13) 


(14) 


At  least  six  unknown  coefficients  are  required  to  be  identified.  By  substituting  the  corner 
frequencies  into  eqn  (6),  we  have  the  identified  transfer  function 

1.000029  -f-  0.206648s  +  0.05966s2 
f+  0.2549245  +b.088779s}  +  0.008985s  ,,5) 

The  comparison  of  the  frequency-response  data  of  eqns  (12)  and  (15)  is  shown  in  Table  2.  These 
results  are  also  satisfactory. 


5.  CONCLUSION 

A  simple  method  has  been  presented  for  fitting  a  transfer  function  from  experimental 
frequency-response  data.  A  logarithmic  amplitude-frequency  curve  is  first  plotted  from  the 
available  frequency-response  data,  then  it  is  smoothed  and  approximated  by  piecewise  seg¬ 
ments  with  integer  or  fractional  slopes.  As  a  result,  the  most  important  interpolation  data  and 
the  order  of  the  transfer  function  may  be  obtained  from  the  irrational  transfer  function.  When 
the  slope  at  two  consecutive  low  frequencies,  <ui  and  <u2,  is 


x(slope)  = 


iro'wi)U  - 


r(/wr)U 


20  log 


(In  other  words  there  exists  x  poles  at  the  origin.),  then  the  available  frequency-response  data 
should  be  multiplied  by  (/»)*  so  that  eqn  (6)  may  be  applied.  The  method  presented  in  this 
paper  is  useful  for  digital  computation  and  provides  an  additional  tool  for  system  identification. 
A  computer  program,  based  on  the  approach  discussed,  has  been  written  in  the  appendix. 
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APPENDIX 

This  program  is  used  to  fit  a  transfer  function  using  frequency-response  data.  The  details  to  prepare  the  input  cards  can 
be  summarized  as  follows: 

The  first  data  card: 

NDT— number  of  available  frequency-response  data 

NP— number  of  different  transfer  function  structures  to  be  identified. 

The  second  data  card: 

XWhj  =  I  to  NDT —a  vector  of  the  frequency  values  at  which  there  is  available  data 
The  third  data  card: 

XRf.j  =  I  to  NDT— a  vector  of  the  values  of  the  real  parts  of  the  available  data  at  XW,. 

The  fourth  data  card: 

X/j,  j  =  I  to  NDT— a  vector  of  the  values  of  the  imaginary  parts  of  the  available  data  at  XW,. 

The  fifth  data  card: 

m— The  number  of  the  unknown  constants  in  the  numerator  polynomial  of  the  transfer  function  to  be  identified. 
n— The  number  of  the  unknown  constants  in  the  denominator  polynomial  of  the  transfer  function  to  be  identified. 

The  sixth  data  card: 

NDj.j  m  1  to  NAf— A  subscript  number  is  assigned  to  each  set  of  frequency-response  data.  ND,  is  the  vector  of  those 
subscript  numbers  which  point  to  the  frequency-response  data  set  to  be  used  to  identify  the 
transfer  function.  NM  =  n+m. 

The  numerical  example  in  Example  I  is  used  to  illustrate  the  procedure.  For  the  given  system,  13  (i.e.  NDT  =  13) 
frequency-response  data  are  available  in  Table  I.  Various  combinations  of  the  structures  of  the  numerator  and 
denominator  polynomials  may  result  in  various  kinds  of  transfer  functions,  if  we  are  interested  in  only  one  (i.e.  NP  =  I) 
structure  of  the  transfer  function,  then 


T(s)  = 


po+PiS 

T+  q,s  +  q7s* 


(16) 


The  data  on  the  first  data  card  are  NDT  =  13  and  NP  =  I  The  values  of  the  frequencies,  real  parts,  and  imaginary  parts 
of  the  available  data  are  given  in  Table  I.  Therefore,  the  data  on  the  subsequent  input  cards  are 


XW,  =  0.I. 
XR,=  1.0000, 
XI,  =  -0.1130. 


XW-  =  0.2 . 

XR;  -  1.0000, . 
X/;  —  —0.2200, 


XW„=IOO 
XR„  =  -0.I130 
X/„  =  0.9930 


The  data  on  the  next  card  is  the  number  of  the  unknowns  in  the  numerator  and  denominator  polynomials  in  eqn  (16): 


m  =  2  and  n  =  2. 


The  corner  frequencies  (the  most  important  data)  occur  at  XWj  =  0.5.  XW6  =  2,  XW,  =  10,  and  XWM  =  40;  therefore, 
the  values  of  the  selected  subscript  numbers  (i.e..  ND )  are  ND,  =  3,  ND:  =  6,  ND,  =  9,  and  ND,  =11.  These  data  appear 
on  the  last  data  card. 

The  output  of  this  program  is  p0  =  0.99628.  pi  =  -0.991402,  q,  =  0.10053  and  q2  =  0.010072.  Also,  the  real  parts, 
imaginary  parts,  magnitudes,  and  phase  angles  at  available  frequencies  of  the  identified  transfer  function  in  eqn  (16)  are 
calculated  and  printed  for  comparison  with  the  given  data. 

A  listing  of  the  computer  program  is  as  follows: 


C  A  PROGRAM  TO  FIT  TRANSFER  FUNCTION  USING  FREQUENCY-RESPONSE  DATA. 
C 

OnUHt  h  PRECISION  »<5O),XR(S0),XI(50>,XRtl (  SO )  .  XRI 2 ( 50 )  , It  C  10 ) , 

IX (JO, 30), G( to, JO) ,OETN,H< I0),XS,X  ,X»C50) 

D1»Fwsiom  KOI  tt> ) 

complex  cx,rxx,cxxx,C(V>)  ,cx t 

1  ■ » 0  0  REAP  (5,50tl  vOT.NP 


bill 

El)R*AT 

t  IMS) 

RRlTt 

(6,601)  NOT, 

NP 

601 

FORMAT 

C(2X,16f5)> 

REALMS 

,S02HX»(J)  , 

J*1  , 

NOT) 

50  2 

FOR  VAT 

(4F70.8) 

REAIUS 

,  S02  )  (XR(.I) 

,  J*1 

,'UT) 

R  K  A  O  (  5 

,502)  (XI(J> 

,J*1 

.NOT) 

DO  10 

J*l ,«UT 

10 

RRITE 

(6,602)  J.X 

R(J) 

, XR ( J ) 

602 

FflHRAT 

(SX, IS.5F20 

.8) 

Ilf)  <tn 

FR=1 , NP 

HEAR  ( 

5,50  1  )  *,N 

*fR  ITE 

(6, out)  N,N 
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NMsNtM 

READ  (5, son  (ND(J),J»1.NH> 

UDITi:  (6,601)  (N[>(J),J»1.NN) 

00  BO  J=1,N* 

JJxNOC J ) 

«(J)=X«(JJ) 

mRITF  (6,602)  JJ , XN ( J J ) , XR ( J J) , X I ( JJ  ) 
XRU(J)»XR(JJ)-XI(JJ) 

80  XRI2(J)*XR(JJ)*XI(JJ) 

00  20  Ksl,N« 

B(K)»XRI2(K) 

A(K.I)«I. 

IF  (N.EQ.l)  GO  TO  21 

LX*  1 

XSel. 

00  30  J*2 ,  M 
A(K,J)*XS*o(K)«*(J-l) 

LK=LK*1 

IF ( LK . EQ . 2  )  XS«(-1.)*XS 
30  1F(LK.E0.2)  LX»0 

21  CONTINUE 
LK«t 
XS=-1  . 

JBSM+l 

DU  40  JeJM.NN 
IF  (I.K.EO.l)  X*XRI1(K) 

IF  (LK.F0.2)  X»XRI2(K) 

(F  (LK.E0.2)  X8«(»1.)*XS 
A(K  ..l)=XS*X*a(K)**(  J-JN  +  U 
LKsLR+l 

40  IF(LK.GT.2)  LK*1 

20  CONTINUE 

CALL  INVER  (A,NM,G,0,DETN,B,H) 

„R1TE  (b,bOI)  n,N,(H(.I),J»1,WN> 
oOI  FORMAT  (//2X,215/(2X.5E20.8)//> 

DO  SO  K-l.NDT 
XXSaH(l) 

CX*CMPLX(XXS,0.) 

XX*XR(R) 

IF(M.EQ.l)  GO  TO  bl 
DO  60  J*2,« 

cxy=cnflx(o, ,xx) 

CXX=H( J)*CXY**J1 
60  CXsCX+CXX 
h|  CXXX*CMPLX(1.,0.) 

DO  70  JsJM.NN 
CXY«C«Pt,X<0.  ,XX) 

CXX*H(J)*CXY»»JJN  . . .  .  .. 

70*  CXXX«CXXX*CXX 
C(K)=CX/CXXX 

SO  NRITE  (6,604)  K,XN(K),C(K) 

604  FOHHAT  (2X,IS,F20.8,tOX,F20.8,5X,F20.H) 

DO  81  J*l,NDl 

XCT*XR(U) 

YCT*XI(J) 

CX  Y*('MPLX  (  XCT ,  YCT  ) 

XNlsCAHStCXY ) 

XT1*ATAN2(YCT,XCT)  *67.2958 
CX»C(J) 

Xm2*CARS(CX) 

XXCTnREAL ( CX ) 

YYCTsAlMAGICX ) 

XT2*ATAN2(YYC'T,XXCT)  *67,2958 

NRITE  (6,605)  J,XN(J),XH(J),Xl(J)  ,XM1,XT1 

81  NKITE  (6,606)  U  .  XN  (  J  )  ,C  U  )  ,XM2,XT2 

605  FORMAT! /2X, I5,F20.8,2X,F20.8,2X,F20.8,2X,E20.8,2X,F20,8) 

606  FORMAT! 2 X, IS »  F20.B,2X,F20,8,2X,F20.8,2X,C20,8,2X,F20,8/) 
90  CONTINUE 

GO  TO  1000 
END 

SUBROUTINE  INVER  ( A , N ,8 , N ,OET , XC , XO) 

DOUBLE  PRECISION  A (  30 , 30 ) ,B( 30 , 10 ) , IPVOTC 30) , INDEX (30 ,23 , 
1PIV0T( 30),XC(30),XD(30) ,DET,T  ,5 
EOUI VALENCE  ( I  HON , JROa ) , ( ICOL, JCUL) 

500  FORMAT  (12) 

501  FORMAT  (4F20.6) 
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601  FORNAT  ( /// ( 2X , 8F1 5 .6 ) > 

M»0 

57  OET*l. 

1)0  17  J=1,N 

17  1 PVOT ( J ) so 
0"  I  16  1*1  ,N 
T  =  0. 

DO  9  J»1,N 

IF(1PVUT(J)-1)  1 3,9, 13 
1  1  DO  23  H  =  1 ,N 

1F(IPV0T(K)-1 )  41, 23, hi 
43  IF  (l>ABS< T)-DABS( A ( J ,K ) ) )  83,23,23 
83  I ROW* J 
1C0L=K 
T*A(J,K) 

23  CONTINUE 
9  CONTINUE 

1PVOTC I COL ) * IPVOT (ICOL)Fl 
1F(IR0W-IC01.)  73,109,13 
73  DETs-OET 
00  12  L*  1  ,  N 
T  =  A(TRU»,L,) 

A  (  IHOw,l,)sA  (  ICUL.L) 

12  A ( 1CUL ,  L  1  =T 
IF(8)  104,109,33 

13  DO  2  L*l,* 

T»8(  IRO» ,  L  ) 

8(1  ROW  ,  L  )  =  R  ( 1CUL  ,  1. 1 

2  R( ICOL , L ) =T 

1U9  INDFX(I,l)=(ROW 
I  NOEX  (  T  ,21  =  ICOC 
Pt  VOT(  I  )  sA  (  ICOL,  ICOI. ) 

DFT=OFT»PIVOT( 1 ) 

A(  icnr.,  i  cm,  )*i. 

on  205  (,*)  ,N 

205  A(ICni,,l,)*A(IC0I,,L)/PIVOT(I) 

IF(N)  347,347,66 
66  DO  52  L*1,H 

52  8(1  COL  ,  1.3*8  <  I  COL  ,  L  )  /  PI  VOT(  I  ) 

347  DO  134  I.Isl.N 

IF  (I.  I -ICOI.  1  21,134,21 
21  T  =  »U.1,1C0L) 

A(LI,IC0L)*0. 

DO  89  L*  1 ,  N 

89  A((,I,l,)*AUI,l.)-A(lCOl,.L)*T 
IF(h)  134,134,18 

18  00  68  [,=  1,« 

68  B(I,I,L)  =  H(LI,L)-H(1C0I,,L)*T 

134  CONTINUE 

135  CONTINUE 
222  DO  3  1*1, N 

L  =  N-I»1  .....  ...  -  .  • 

1  F (  INDEX (L,  1  )-INDt'X(L,2)  )  19,3,19 

14  JRO.sINOEXU,,  1  ) 

jcol* index ( l , 2 i 

00  549  8*1, N 
T*A(F,JR0W) 

A(K,J8U«)=A(K,JC0I) 

A(8,JCOI.)sT 
549  CONTINUE. 

3  CONTINUE 

00  40  K* 1 , N 
40  CONTtNUF. 

DO  20  K*1,N 
S*0. 

DO  30  J*I,N 
30  S«Se A(K,J)*XC(J) 

20  XD(K)*S 

•RITE  (6,601)  (X0(8),F*l,N) 

81  CONTINUE 
RETURN 
END 
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Solution  of  state-space  equations  via  block-pulse  functions 

L.  S.  SHIEHf,  C.  K.  YEUNGf  and  B.  C.  McINNISf 

A  recursive  algorithm  is  developed  for  the  piecewise -constant  solution  of  dynamic 
equations  via  block-pulse  functions  where  j=  1,  2,  m.  For  l^j^m  (where 

j  anti  m  are  integers)  and  final  time  T,  each  block-pulse  function  fa{t)  *8  defined 
by  fa(t)=  1  for  (j  —  \)TJm^t  <  jTJm  and  =  0  otherwise.  Compared  with  Walsh 
function  approaches,  the  proposed  method  is  simpler  to  compute,  is  more  suitable 
for  computer  programming,  and  provides  the  same  accuracy.  Also,  a  discrete-time 
solution  is  derived  for  a  zero- input  state  equation. 

1.  Introduction 

Consider  a  linear  time-invariant  system  described  by  the  state  equation 

x(t)  =  Ax(t)  +  Bu(t)  (la) 

and  an  initial  vector 

*(0)  =  *o  0  b) 

where  A  is  an  n  x  n  system  matrix,  B  is  an  n  x  r  constant  matrix,  x(t)  is  a  state 
vector  of  n  components,  x(l)  is  a  rate  vector,  and  u(t)  is  an  r-component  input 

I 

vector.  It  is  often  difficult  to  evaluate  the  integration  J  x(t)  dl,  which 

0 

is  the  solution  x(t)  in  (1),  by  a  numerical  method  (Carnahan  et  al.  1969).  One 
approach  is  to  find  a  set  of  orthogonal  functions  for  the  approximate 
solution  as  follows  : 

t  i 

x(t)=x(0)+  J  x(t)dt*P  j  </'(t)dtzPQf(t)=WMt)  (2) 

o  o 

where  P,  Q  and  W  are  n  x  m,  m  x  m  and  nxm  weighting  matrices,  respectively, 
and  is  an  mx  1  vector  with  m  orthogonal  functions  ^(t),  which  are  both 
suitable  for  approximation  of  x(t)  and  easy  to  integrate  numerically. 
Corrington  (1973),  Chen  and  Hsiao  (1975),  and  Rao  and  Sivakumar  (1975) 
chose  Walsh  functions  as  the  tp ((t)  for  the  approximate  solution  in  (2) 
and  reported  that  their  piecewise-constant  solution  gives  a  satisfactory 
result.  However,  their  computational  methods  (Chen  and  Hsiao  1975,  Rao 
and  Sivakumar  1975)  either  required  the  inversion  of  a  large  matrix  or  the 
inversion  of  many  small  matrices.  This  results  in  computing  time  and  storage 
being  wasted,  and  the  truncation  and  round-off  errors  might  be  seriously 
accumulated.  Recently,  Chen  et  al.  (1976)  and  Gopalsami  and  Deekshatulu 
(1976)  introduced  a  set  of  ‘  block-pulse  functions  ’  for  the  solutions  of  distri¬ 
buted  systems  and  identification  problems.  They  pointed  out  that  there  is 
a  one-to-one  relationship  between  Walsh  functions  and  block-pulse  functions. 
For  I  where  j  and  m  are  integers,  the  block-pulse  function  is 

» 
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re-defined  and  extended  in  the  interval  0  $  t  <  T  (rather  than  in  the  interval 
0</<  1  as  in  Chen  et  al.  (1976),  and  Gopalsami  and  Deekshatulu  (1976))  and 
by 

(1  for  (j  —  1  )T/m < t<j Tjm 

(3  a) 

0  otherwise 

T  is  the  final  time,  and  m  is  the  number  of  subintervals  between  <  =  0  and 
t  —  T  as  well  as  the  number  of  block-pulse  functions  to  be  used.  When  m 
block-pulse  functions  are  used  to  approximate  the  integration  of  the  original 
block-pulse  functions,  we  have 


I  rp  rp 

f  4,(1)  dl^-  H4>(t)  =  - 

o  m  m 


■j  i  i  •  r 

'  M) ' 

0  J  1  •  1 

4>S) 

_o  o  o  • 

.4>Jt). 

(3  6) 


where  <£(<)  is  an  m  x  1  vector  with  m  block-pulse  functions.  The  constant 
matrix  (Tjm)H,  with  the  dimensions  mxm,  is  the  operational  matrix  (Chen 
et  al.  1976,  Gopalsami  and  Deekshatulu  1976)  for  the  block-pulse  functions. 
Sannuti  (1976)  discussed  the  properties  of  the  <^(f)  and  proposed  a  method 
for  the  solutions  of  linear  and  non-linear  problems.  From  (3  b)  we  observe 
that  the  matrix  H  is  an  upper  triangular  matrix  that  consists  of  diagonal 
elements  being  £  and  the  other  elements  being  1.  By  taking  advantage  of 
this  peculiar  arrangement  of  H  and  by  choosing  the  block-pulse  functions 
4>At)  as  the  ^r(l)  in  (2),  an  alternative  method  is  proposed  in  this  paper 
to  derive  an  effective  algorithm  for  the  piecewise-constant  solution  of 
the  state  equation  in  (1).  The  computation  in  our  algorithm  involves  the 
inversion  of  only  one  matrix  that  has  the  same  size  as  the  system  matrix. 
Compared  with  Walsh  function  approaches  (Corrington  1973,  Chen  and  Hsiao 
1975,  Rao  and  Sivakumar  1975)  the  proposed  method  is  simpler  to  compute, 
is  more  suitable  for  computer  programming,  and  provides  the  same  accuracy. 


2.  Main  result 

Let  X;(t)  be  the  ith  component  of  the  state  vector  x(t)  that  is  the  solution 
of  the  state  equation  in  (1).  The  x,(l)  can  be  expressed  approximately  as 

tn 

£  C,  )4>j(t),  where  m  is  a  large  finite  number,  4>j(t)  are  block-pulse  functions, 

j  =  \ 

and  Gtj  are  weighting  constants  to  be  determined.  The  state  vector  x(t)  can 
also  be  approximated  as 

x(t)=zC4>(t)  (4  a) 

where 

C  =  [CvCt . C„]  (4  6) 

and 

. *■<*)]'  (^) 

The  prime  designates  the  transpose,  and  the  nxm  matrix  C  consists  of  m 
column  vectors  Cj  to  be  determined.  Our  goal  is  to  develop  an  effective 
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algorithm  to  determine  Cj  for  every  j  so  that  the  piecewise-constant  solution 
in  (4  a)  can  be  obtained. 

We  will  now  derive  the  recursive  algorithm.  Let  the  rate  vector  i(t)  in 
(1)  be  approximated  as 

x(t)^D<j>(t)  (5  a) 

by  using  m  block-pulse  functions,  where 

D  =  [dj,  d2»  •••>  ]  (5  b) 

The  D  is  an  n  x  m  constant  matrix  with  m  column  vectors  dj  of  size  nx  1  to 
be  determined.  Integrating  (5  a)  and  using  the  results  of  (3  b)  and  (4  a) 
yields 


where 


x(t)^D  C  <f>(t)dt  +  x( 0)s  -  DH  +  G  <f>(t)  =  Cd>(t) 

o  Lm  J 

G  =  [#(0),  x(0) . a;(0)]  =  [gx,  gt . grm] 


C =— DH  +  G  =  [CV  C2,  Cm]  (6  c) 

771 

The  gt  in  (6  6)  is  the  initial  vector  x(0),  and  the  constant  matrix  (Tjm)H  is 
shown  in  (3  b).  The  accuracy  of  an  approximate  solution  in  (6  a)  depends 
on  the  number  of  block-pulse  functions  and  the  time  interval  Tjm  used.  The 
rx  1  input  vector  u(t)  in  (1)  can  also  be  approximated  as 

u(t)^L<f>(t)  (7  a) 

using  m  block-pulse  functions,  where 

L  =  [LvL2 . Lm ]  (7  b) 

The  rxm  matrix  L  consists  of  m  column  vectors  L}  to  be  determined.  Bv 
applying  the  orthogonal  property  of  the  block-pulse  functions  to  (7  a),  we 
have 

yyi  jT  I  tn 

Lj  =  r  f  w(t)dtS^[u(jTjm)  +  u((j-l)Tlm)]  (7  c) 

•»  (j-l)T/m 

equals  average  vqlue  of  u(l)  over  the  interval  (j  —  1  )T/m  $  t  ^  jT/rn.  The 
accuracy  of  the  approximation  in  (7  c)  depends  on  the  time  interval  Tjm 
used.  Substituting  (5  a),  (6  a)  and  (7  a)  into  (1  a)  yields 

D-  —  ADH  +  AG  +  BL  =  —  ADH  +  K  (8a) 

m  m 

where 

K  =  AG  +  BL  =  [kvk . JfcJ  (8  5) 

The  column  vector  k}  is  an  n  x  1  known  vector.  The  unknown  matrix  D  in 
(8  a)  and  (5  a)  can  be  determined  from  the  matrix  equation  (Chen  and  Hsiao 
1975) 


A  ® m  77  J  —  m  J  —  (®  c) 
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or 

.4,  0  <(  ...  0 

-A  A,  0  ...  0 

-A  -A  At  ...  0 

-A  -A  -A  ...  Ai 

where 

(He) 

The  lnm  in  (8  c)  is  an  tint  x  nm  identity  matrix,  and  the  (g>  in  (8  c)  represents 
the  Kronecker  product.  Each  nxn  biock  element  0  in  (8  d)  is  a  null  matrix 
and  /„  in  (8  e)  is  an  n  x  n  identity  matrix.  It  is  known  that,  as  more  orthogonal 
functions  are  used  to  approximate  x(t),  a  better  approximate  solution  is 
obtained.  Therefore,  m  should  be  a  large  number  and  the  matrix 

((mlT)I„ul-A®H‘] 

is  large.  The  direct  inversion  of  such  a  matrix  for  the  solution  of  d,  in  (8  r) 
is  not  an  effective  method  as  far  as  the  computing  time  and  storage  are 
concerned.  However,  from  the  peculiar  formulation  of  the  square  matrix 
in  (8  d),  we  can  derive  an  effective  algorithm  for  solving  dj  instead  of  invert¬ 
ing  the  matrix  directly.  This  effective  algorithm  is  derived  in  the  following 
way.  By  pre-multiplying  each  block  element  on  both  sides  of  (8  d)  by 
and  by  rearranging  the  new  matrix  equation,  we  have  an  alternative  form  of 
(8  d)  as 

[dti  r/?2  6  0  ...  o i r  d,  i  r/?Al 


where 


Equation  (9  a)  can  be  solved  readily  for  dj.  After  obtaining  the  matrices  if, 
and  f?2  and  the  vector  d,  in  (9  6),  we  can  immediately  determine  the  vector 
d2  from  the  first  equation  in  (9  a).  Then  we  can  substitute  d2  into  the  second 
equation  and  solve  for  d3,  etc.  Note  that  the  m  can  be  chosen  so  that 
((mlT)I„-(\)A)~1  exists. 


i?2  i?2  0  ...  0  d2  w  if)^3 

ff.j  i?2  R2  ...  0  d3  1  R1ki 


R2  i?2  if 2  •••  if2 J  |_dm_,_ 


R^Ar'-^I^-lAy 
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The  general  algorithm  is 


m 

®1  =  ~ji  -®1^1 


where 


7/1  ^ 

dj  =  ^2  X  +  R\kj  =  dj_l  +  y  R\(kj  —  kj-i) 

for  j=2,  3,  ....  »i 


*i=(9 


Rt  =  A1~l  A=RXA  J 

Consequently  from  (6  c)  and  (10)  we  have  the  required  column  vectors  Cjt  or 

°i='Ldi+gi 

i-i  T  (n) 

Gi  =  ~  t  di  +  2^di  +  9i  =  Ci-i  +  2^,  (di-i+di)  for  J=2>  3>  ->  m 

Substituting  (11)  into  (4)  yields  the  required  piecewise-constant  solution  of 
the  state  equation  in  (1).  Note  that  the  <f>j(t)  differs  from  zero  only  in  the 
interval  (j  —  l)Tjm^t <  jTjm  ;  therefore,  the  jth  column  vector  C;  is  the 
required  piecewise-constant  solution  in  that  interval.  Another  advantage  of 
the  proposed  method  is  that  Cj  involves  only  the  vectors  d„  k,  and  git  for 
i  =  1  ...  j,  whereas  the  Walsh-function  approaches  (Corrington  1973,  Chen 
and  Hsiao  1975,  Rao  and  Sivakumar  1975)  require  a  whole  matrix  M'  and  a 
whole  vector  \fi(t)  in  (2). 

If  «(<)  =  0  in  (1),  (10)  and  (11)  can  be  expressed  by  a  set  of  difference 
equations 

d(\)  =  ~R^c(  0)  (12  a) 

d(j  +  \)  =  (In+Rt)d(j)  for  j  =  l,  2,  ....  m-  1  (12  6) 

and 

c(\)  =  \(2ln  +  R2)x(0)  (13  a) 

C(j  +  1)  =c(j)  +  ^  (21  n  +  RJdV)  for  j  =  1,  2,  ...,  w—  1  (136) 

The  solution  of  (12)  is 


d(j)  =  (/„  +  R2y-id(  1 )  =  y  (/„  +  RtY-'Brf  0) 

Substituting  (14)  into  (13  6)  yields 

C(1 )  =  J(2/„  +  R2)x(0) 

c{j  +  1)  =  e(  j)  +  mn  +  R2)(In  +  R2)>-'R*c{  0) 
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The  solution  of  (15)  is 

c(j  -M)  =  c(l)  +  \(2In  +  R2)  '£  (7n  +  i?2)^^(0) 

i  =  0 


=  ^(2^n+  #2)  + 

Since  the  trapezoidal  rule  (as  shown  in  (7))  is  used  as  a  base  for  the  numerical 
integration,  or 


I'  (In+RJX* 


x(0)  for  j  =  1,2,...  (16) 


c(j  +  l)  = 


**(j  +  l)  +  **(j) 
2 


(17  a) 


where  x*(j )  is  the  discrete-time  solution,  therefore 

x*(j+i)=  —x*(j)  +  2c(j+  1)  (17  fc) 

Substituting  (16)  into  (17  6)  we  have  the  required  discrete-time  equations 
z*(0)  =  x(0)  } 


x*(l)  =  (In  +  R2)x(0) 

x*(j+l)=-x*(j)  +  (2In  +  R2)[ln  +  (/„  +  /f2)^2Jz(0) 


}  (18a) 


The  solution  of  (18  a)  is 

x*U)  =  (/„  +  Rt)>x( 0)  for  j  =  0,  1 ,  2,  .. .  (18  b) 

(m  \-»  T 

where  R2  =  [—In  —  \A\  A,  T  =  the  final  time,  and  the  sampling  period  =—. 

Equation  (18  6)  can  be  further  analysed  as 

x*(j)  =  Un  +  **]'*«>)  =  <t>*(j)x(0)  for  j=0,  1,2,  ...  (19a) 

where 


<t>*(j)  =  the  transition  matrix  of  a  discrete-time  system 


~[^n+  RiV 

=  [In  +  (In-lAAT)~lAATy  for  J=0,  1,2 . and  A 5T  =  - 

m 

=[(/„-MAT)-M4+MAT)V 

=  [/„+  4AT+i(AAT)8  +  ^  (AAT)»  +  i  (AAT)*  + 


=  [ln  +  AAT  +  }(AAT)°  +  £J-i(AAT)<J 

The  exact  solution  of  (1)  (with  u(t)  =  0)  is 

x(/)=exp  (^4<)x(0)  =  «I>(<)*(0) 


(19  6) 


(20  a) 
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where 

<!>(()  =exp  (At)  =  the  transition  matrix  of  a  continuous-time  system 
=  [exp  (-4  A  ST)]-#  for  j  =0,  i,  2,  3,  and  t  =jAT 

=  ^In  +  AAT  +  j[(AAT)‘i  +  ^  (/lAT)3-^  (AAT)*+... 

In  +  AAT  +  l(AAT)2+  £  i(-4AT)<T  (206) 

i  =  3  *'  J 

Comparing  with  <t>(t)  we  observe  that  the  first  three  terms  of  (19  6) 

are  equal  to  those  of  (20  6),  while  other  terms  differ  in  weighting  factors 
l/2i_1  in  (196)  and  1/i!  =  l/i(i  —  l)(t  —  2) ...  1  in  (20  6).  Therefore  <f>*(^')  is  a 
good  approximation  of  <D(<)  if  AT  is  small.  Also,  we  observe  that  <I>*(j)  is 
a  finite  matrix,  while  <P(t)  is  an  infinite  series  of  matrices,  therefore  it  is  more 
convenient  to  evaluate  than  $>(<). 

It  is  believed  that  the  derivation  of  the  approximation  of  <t>(t)  in  (20  6)  by 
0>*(j)  in  (19  6)  is  new.  When  u(t)  #  0,  the  approximate  discrete-time  solution 
x*(t)  of  x(t)  in  (1)  can  be  obtained  from  (11)  and  (17  6). 


3.  An  illustrative  example 

Consider  the  dynamic  equation 

x(t)  =  Ax(t)  +  Bu(t) 
x(0)=x0 

where 


~1  2' 

'2  O' 

"l” 

A  = 

,  B  = 

,  *(0)  = 

3  -4 

1  1 

1 

and 

u(t)  —  unit-step  functions 


(21) 


♦  l 

1  /li  i 

/♦l  dt 
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.  ...t 

°  I/*. 

1*2 
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1  /L 

/*2  dt 
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ru 

r  o 
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Figure  1.  The  block-pulse  functions  and  their  integrations. 
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The  piecewise-constant  solution  of  the  state  equation  is 

x(t)sC<f>(t)  (22) 

The  block-pulse  functions  <f>}(t)  and  the  integration  of  the  are  shown  in 
Fig.  1.  The  C  is  an  unknown  matrix  to  be  determined.  The  steps  to  deter¬ 
mine  C  can  be  listed  as  follows  : 

Step  l 

Choose  T  =  1  s  and  in  =  4.  This  means  that  four  block-pulse  functions 

<f>j(t),  j  =  1 . 4,  are  used  in  the  interval  $  1,  and  the  sampling  period  = 

Tjm  =  0-25s. 
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Step  5 

Evaluate  the  required  C  in  (11), 

C  =  \CV  C2,  C3,  C\\ 

where 


c,~d,+s<- 

'1-7949' 

1-2820 

T 

'3-8773' 

C1,  —  Ci  +  —  (d2  +  dt)  — 

2-1792 

O' 

"7-3038" 

C3  =  C2  4-  —  (d3  +  d2)  = 

3-8647 

'T 

'13-0125 

^4  -  C3  +  —  (d4  +  d3)  - 

6-6936 

The  required  piecewise-constant  solution  in  (21)  is 

xt(t)  2  1-7949^(0  +  3-8773^2(0  +  7-3038^3(0  +  13-0125^(<) 
xt(t)  S  1-2820^(0  +  2-1792^,(0  +  3-8547<£3(0  +  6-6936 ^(t) 


Figure  2.  The  exact  solutions  and  the  approximated  solutions. 
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The  exact  solution  of  (21)  is 


*1  (t)  =  ¥  exp  (2t)  -  &  exp  ( -  5t )  - g 
xt(t)  =  |  exp  (2<)  +  &  exp  ( -  5t)  -  f 

The  response  curves  of  the  exact  solution  and  the  approximated  solution  are 
shown  in  Fig.  2.  The  approximate  discrete-time  solution  x*(t)  of  x (t)  in  (21) 
can  be  obtained  from  the  C  in  (22)  and  (17  b). 

If  w(/)  =  0  in  (21),  the  exact  solution  of  (21)  is 


xt(t)  =  |  exp  (2<)  —  f  exp  ( -  51) 
x2(t)  =  |  exp  (21)  -)-  $  exp  ( —  51) 


(23) 


From  (23)  and  (18)  we  can  evaluate  the  exact  solution  x(t)  and  the  approxi¬ 
mated  solution  x*(t)  at  samples  j  =  1,  2,  3,  4,  and  sampling  period  —  T/m  — 
0-25.  The  results  are  tabulated  as  follows  : 


j 

t 

(0 

Zi*(f) 

*z(l) 

***(<) 

0 

0 

1 

1 

1 

1 

I 

0-25 

1-843 

1-872 

1-051 

2 

0-50 

3-095 

3-167 

1-589 

1-610 

3 

0-75 

5-119 

5-289 

2-571 

2-651 

4 

100 

8-444 

8-818 

4-225 

4-411 

It  is  interesting  to  observe  that  the  solution  obtained  by  the  four-point 
approximation  is  quite  satisfactory. 
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Abstract— A  method  is  given  for  optimally  fitting  parameter  matrices  of  state  equations  from  the  real  and 
imaginary  parts  of  noise  free  frequency-response  data  of  a  multi-input,  multi-output,  linear  dynamic  system. 
It  is  assumed  that  all  state  variables  are  accessible  for  measurement.  The  obtained  data  may  contain 
measurement  errors. 


I.  INTRODUCTION 

Several  authors  [1-3]  have  considered  the  application  of  frequency  response  concepts  for 
identification  of  dynamic  systems.  The  problems  of  predicting  parametric  error  from  frequency 
response  measurements  have  also  been  investigated  [3-5],  A  method  is  presented  here  to 
determine  the  best  estimate,  in  least  mean  square  sense,  of  the  parameter  matrices  of  the 
multi-input,  multi-output,  linear,  time-invariant  dynamic  system  equations  if  all  the  state 
variables  are  accessible  for  measurement.  The  obtained  data  are  noise  free  and  contain 
measurement  errors. 


2.  DERIVATION 

The  state  equations  of  an  asymptotically  stable,  completely  controllable  and  observable, 
linear  time-invariant  system  are  given  by: 

£  =  AX  +  BU  (1) 

Y  =  Cx  (la) 

*(0)  =  [0J  (lb) 

where  A  is  a  constant  n  x  n  system  matrix,  X  is  an  n  x  1  state  vector,  B  is  a  constant  n  x  r  input 
matrix,  C  is  a  constant  mxn  output  matrix,  U  is  an  r x  1  input  vector,  and  Y  is  an  m x l 
output  vector.  Let  us  define. 


B  =  [fc„...,fir] 


where  Si  is  an  n  x  1  column  vector  and  C,T  is  an  n  x  1  row  vector. 


(2) 

(2a) 


(2b) 
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The  Laplace  transformation  of  eqns  (1)  and  (la)  yields, 

(s/  -  A)X(s)  =  BU  (s)  (3) 

and 

?(s)  =  £*(s).  (3a) 

Successive  choice  of  each  of  the  scalars  L/t(s)  in  eqn  (2b)  as  an  input  while  the  remaining 
scalar  components  of  C(s)  are  zero  yields  the  following  set  of  transfer  functions  from  each  of 
the  scalar  inputs  to  the  state  variables. 


(sf-A)t(s)  =  bt 


(4) 


where 


f,(s)  =  £P^&(r)  and  e=\ . r.  (4a) 

If  the  input  functions  of  U(s)  are  sinusoidal  functions  with  varying  frequencies  ««*,  we 
obtain  the  corresponding  frequency  response  data  f,(ja»*)  as  follows: 

fe<jo>k)  =  Pt(m)+ i4'(a>k)<  e=  1 . r  (5) 

where  P,(<o t)  and  4«(<m)  are  vectors  of  the  real  and  the  imaginary  parts  of  fk(jtok). 

Multiplying  the  steady  state  portion  of  eqn  (4a)  by  a  normalization  constant  M,  (i.e.  the 
magnitude  of  a  sinusoidal  input  function)  we  have 

X,(ja>k )  =  )  =  MJPAvk )  +  /M4,(<m ),  e  =  I . r  (5a) 

and 

*(M)=  2  *«(/<*>  <5b> 

Substituting  s  =  juk  and  eqns  (5)  and  (3a)  into  eqns  (4)  and  (3a)  yields 

l jo>kf  -  -4][P,(«i»»)  +  j4'M)  =  b,  (6) 

Y,(j<ak)  =  C\Mf,(wk)  +jM'4,(u>k)]  ~  &(<m)  +  /M«m)  (6a) 

and 

Yih>k)  =  2  fcO**)  «>b> 

e-l 

where  £,(<m)  and  Me*)  are  vectors  of  the  real  and  imaginary  parts  of  Y,(jtok).  After  we 
separate  the  real  and  imaginary  parts  of  eqns  (6)  and  (6a)  and  equate  the  respective  real  and 
imaginary  parts,  we  have 


A4Ao»k)  =  <okP,(a>k) 

(7) 

ApA&k)  +  cm 4«(«m)  = 

(7a) 

#r(«m)  =  CM'PAak) 

(7b) 

(7c) 

The  parameter  matrices  A,  6,  and  C  can  be  obtained  as  follows: 


t 
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A  =  [oil P,(«i).  <tl2p,((i>2),  ....  W.A(«.)J[4.(«»|).  4»("2) . 4«(««)]"'  (8) 

b,  =  -(AP,(o»*)  +  <»ir4«(  <*)*)]  (8a) 

and 

C  =  [4,(011),  4,(o»j),  ....  Ot|),  A^A(<U2),  •  •  ■  >  Af^,(o»„)]  1 

=  [Mai,),  Mtii]) . Moi.)][M4,(a>i),  A#,4,(a>2) . Ma(»«)] 

=  [2  8.(0.,), ....  2  8,(oi.)][S  M^.(oi,). ...,2  M^.(o,»)]  ' 

=  [2  M«i),  •  •  • ,  2  M*>» >][2,  H4«(»i) . 2(  M4«(»«) j  •  (8b> 

The  data  in  eqns  (8H8b)  can  be  chosen  so  that  the  matrix  inversions  exist. 


3.  EVALUATION  OF  OPTIMAL  PARAMETER  MATRICS 
If  the  frequency  response  data  are  noise  free  and  measurement  error  free,  then  there  exist 
unique  parameter  matrices  A ,  B  and  C.  However,  in  practice,  there  exist  measurement  errors 
even  if  the  system  is  noise  free.  As  a  result,  the  evaluated  parameter  matrices  have  inaccuracies 
due  to  the  errors.  In  this  paper,  optimal  parameter  matrices  are  evaluated  from  the  measure¬ 
ment  error  contaminated  data. 

Consider  i  sets  of  parameter  matrices  A,  b,  and  C,  which  are  defined  as  A„  b„  and  C„  and 
which  are  evaluated  from  /  sets  of  data  using  one  control  input  or  r  control  inputs.  If  many  sets 
of  experimental  frequency  response  data  can  be  obtained,  then  the  optimal  parameter  matrices 
A,  b,  and  C  can  be  obtained  from  the  matrix-mean  values,  or 

A  =  t2a.  (9) 

A=r2^  <*> 

C  =  }gC.  (9b) 

However,  to  obtain  many  sets  of  frequency  response  data  is  often  not  practical  and  sometimes 
impossible.  The  following  technique  is  proposed  to  obtain  the  optimal  matrices  with  fewer  sets 
of  frequency  response  data.  Suppose  that  the  system  matrices  Af,  i  =  1, . . . ,  r  can  be  evaluated 
by  r  sets  of  frequency  response  data  which  are  obtained  from  the  controllable  system  by  any 
one  input  U,  or  by  r  sets  of  inputs,  then  we  construct  the  following  matrix  equation. 


where 


£A  =  fi 


A,'1' 

’  M 

Af' 

u 

.A,'1. 

and  P  = 

.  4. 

(10) 


(10a) 


in  which  Af 1  are  n  x  n  inverse  matrices  of  A,  obtained  by  the  use  of  eqn  (8)  and  /,  are  n  x  n 
identity  matrices.  The  desired  optimal  matrix [6, 7]  A  which  minimizes  the  sum  of  squares  of 
residuals  §  =  where  A  =  P  -  PA,  is  given  by 


(ID 


By  a  similar  approach  the  optimal  matrices  b,  and  C  can  be  obtained  as  follows: 
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To  obtain  b,  we  construct  the  matrix  equation 


where 


G,= 


Gr 


La-  J 


02) 


(12a) 


in  which  6,r'  are  n  x  n  inverse  matrices  of  A  and  the  elements  at  jth  row  and  fcth  column  in  C,  and 
H,  are: 


G,0,*)  =  M/,1) 

if  j  =  k 

=  0 

if  j*k 

(12b) 

H'(j,k)=b'(j,  1) 

if  j  =  k 

=  0 

if  j*k 

(12c) 

j  =  1 . n,  /  = 

1 . r 

II 

s: 

II 

l,...,r. 

It  is  interesting  to  note  the  fact  that  </,(;,  k)  and  Ht(J,  k)  are  diagonal  which  greatly  reduces  the 
practical  problem  of  calculating  H* 

The  optimal  matrix  H,  is 


H,  =  (G/Gtr‘  G,TF.  (13) 

The  optimal  vector  b,  can  be  obtained  from  eqn  (12c).  To  obtain  the  optimal  row  vector  C/  in 
C  we  use  the  following  matrix  equation: 


where 


AS,  =  F 


'or1' 

Di' 


D.= 


(14) 


(14a) 


Dr' 


in  which  A*1  are  n  x  n  inverse  matrices  of  A  and  the  elements  of  the  /th  row  and  fcth  column 
in  A  and  §  are 


A(/,*)=CJ(1,/') 

if  /  =  k 

=  0 

if  j*k 

(14b) 

SAbk)=  C/(l,/) 

if  j  =  k 

=  0 

if  j*k 

(14c) 

j  =  1 . n,  i  =  1, . . . ,  r,  z  =  1, . . . ,  m,  k  =  1, . . . ,  n.  Here  again  the  structure  of  AO’.  *)  >s  quite 

favorable  for  performing  the  necessary  inversions. 

The  optimal  matrix  St  can  be  obtained  from 


4-(  ATAr'ATfi 


(15) 
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The  optimal  row  vector  C/  can  be  obtained  from  eqn  (14c).  After  obtaining  the  optimal  vectors 
6,  and  C.T  we  have  the  optimal  input  matrix  and  output  matrix,  or 


and 


(16) 


(17) 


4.  CONCLUSION 

A  method  for  the  solution  of  the  difficult  problem  of  identifying  a  multi-input,  multi-output, 
linear  system  from  measurement  error  contaminated  data  has  been  presented.  The  resultant 
parameter  matrices  are  optimal  in  the  least  mean  square  sense.  The  particular  advantage  of  this 
technique  is  the  ability  to  utilize  a  relatively  limited  amount  of  experimental  data  to  obtain  the 
systems  dynamic  equations.  The  identification  process  can  be  easily  performed  using  digital 
computers. 
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APPENDIX 

Illustrative  example.  For  a  known  dynamic  system  described  by  the  following  state  equation, 

k  =  AX  +  BV 

f  =  CX  (18) 

where 

4-fl^  Mil]  -4-Cfl 

the  error  contaminated  frequency  response  of  Table  1  was  obtained. 

Assuming  a  unity  magnitude  for  the  excitation  function  or  M,  and  M„  equal  unity  and  by  following  eqns  (8),  (8a)  and  (8b), 
we  have 


("-1.018087614  -0.967404614 
L  1.963824775  -  3.934809232 

(-0.993028846  -1.015416667 
l  2.005288462  -4.00249999 


] 

] 


(19) 


Table  I.  Frequency  response  data 


<  =  l. 

e  =  2, 

r,(M) 

e=l. 

r,(M> 

<  =  2, 

n(M) 

ft* 

«i(u») 

/W 

M«») 

gj(«6) 

*2<«6) 

0.2 

0.658 

-0.077 

-0.821 

0.104 

0.82 

-0.10 

-0.91 

0.13 

0.326 

-0.055 

-0.158 

0.060 

0.65 

-0.11 

-0.32 

0.12 

2.0 

0.269 

-0.346 

-0.288 

0.442 

0.29 

-0.44 

-0.20 

0.51 

0.038 

-0.192 

0.173 

0.135 

0.08 

-0.38 

0.35 

0.27 
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for  e  »  1,  we  have 

f  1.000275632] 

*'  “  L-0.00I001583J* 

r  1.000208233  ] 

[o.OOl  1 50508  ij 

(19a) 

and  when  e  =  2.  we  have 

■  r-1.003112556] 

''I  0.989774387 1 

-  f- 1.00220058  ] 

0.987559404 J 

(19b) 

When  z  =  1.  we  have 

M  . 

Cii  ■ 

=  (1.011006541.  0.474716861], 

Cl  =  [1.007961095,  0.521923674] 

(19c) 

and  for  z  =  2,  we  obtain 

cl- 

[0.0220130807,  1.949433722], 

Cf2=  [0.0003199367,  2.023653999] 

(I9d) 

Applying  eqns  (II),  (13)  and  (IS)  we  have  the  optimal  parameter  matrices 

*-4-1005 
L  1.983 

-0.991] 

-3.966J 

(20) 

fl-f1000 

B  [0.000 

-1.002] 

0.988  J 

(20a) 

and 

C  =  f1009 
Lo.ooo 

0.496] 

1.985  J' 

(20b) 

Compared  with  eqn  (18).  the  answer  is  quite  satisfactory. 
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An  algebraic  method  to  determine  the  common  divisor, 
poles  and  transmission  zeros  of  matrix  transfer  functions 

L.  S.  SHIEHf,  Y.  J.  WEIf  and  J.  M.  NAVARROJ 


A  purely  algebraic  method  which  uses  the  matrix  Kouth  algorithm  and  its  reverse 
process  of  the  algorithm  is  presented  to  decompose  a  matrix  transfer  function  into  a 
pair  of  right  co-prime  polynomial  matrices  or  left  co-prime  polynomial  matrices. 
The  poles  and  transmission  zeros  of  the  matrix  transfer  function  arc  determined  from 
a  pair  of  relatively  prime  polynomial  matrices.  Also,  the  common  divisor  of  two 
matrix  polynomials  can  be  obtained  by  using  the  matrix  Kouth  algorithm  and  the 
matrix  Routh  array. 


1.  Introduction 

The  properties  and  applications  of  poles  and  transmission  zeros  of  a  multi- 
variable  system  have  been  extensively  studied  in  recent  years  by  many 
researchers  (Desoer  and  Sehulman  1974,  Kwakernaak  and  Sivan  1972, 
Rosenbrock  1970,  Moore  and  Silverman  1972,  Wolovich  1972,  1973,  Davison 
and  Wang  1974,  Francis  and  Wonham  1975,  Sinswat  et  al.  1976,  Kouvaritakis 
and  MacFarlane  1976,  Wang  and  Desoer  1972).  Desoer  and  Sehulman  (1974) 
defined  the  poles  as  real  or  complex  numbers  for  which  the  responses  of  a 
circuit  or  system  to  a  series  of  singular  inputs  are  purely  exponential.  The 
transmission  zeros  are  also  defined  as  real  or  complex  numbers  for  which  the 
transmission  of  some  particular  signals  is  completely  blocked.  The  role  of 
poles  in  the  analysis  and  synthesis  of  circuits  and  systems  is  well  known,  and 
in  recent  years  the  transmission  zeros  are  found  to  be  important  in  many 
aspects  of  feedback  control  theory  (Desoer  and  Sehulman  1974,  Kwakernaak 
and  Sivan  1972,  Rosenbrock  1970,  Moore  and  Silverman  1972,  Wolovich 
1972,  1973,  Davison  and  Wang  1974,  Francis  and  Wonham  1975,  Sinswat  et  al. 
1976,  Kouvaritakis  and  MacFarlane  1976,  Wang  and  Desoer  1972).  Therefore, 
it  is  useful  and  desirable  to  have  an  effective  method  to  determine  the  locations 
of  these  poles  and  transmission  zeros.  Several  methods  are  available  to  locate 
the  positions  of  these  poles  and  zeros  (Kwakernaak  and  Sivan  1972,  Rosenbrock 
1970,  Moore  and  Silverman  1972,  Wolovich  1973,  Davison  and  Wang  1974, 
Francis  and  Wonham  1975,  Sinswat  et  al.  1976,  Kouvaritakis  and  MacFarlane 
1976).  However,  most  of  the  suggested  approaches  (Rosenbrock  1970, 
Moore  and  Silverman  1972,  Wolovich  1973,  Davison  and  Wang  1974,  Francis 
and  Wonham  1975,  Sinswat  et  al.  1976,  Kouvaritakis  and  MacFarlane  1976) 
are  derived  for  the  systems  which  are  represented  by  state  equations  in  the 
time  domain.  The  major  disadvantage  of  most  time-domain  approaches  is 
that  the  computation  may  not  be  very  attractive  if  the  dynamic  systems  are 
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of  high  order.  When  a  given  multivariable  system  is  described  by  a  matrix 
transfer  function  that  might  have  a  high  degree  common  divisor  (the  common 
factor)  of  the  numerator  and  denominator  matrix  polynomials,  the  order 
of  the  corresponding  state  equations  is  in  general  very  high.  Therefore, 
most  time-domain  approaches  may  be  difficult  to  apply.  In  this  paper,  a 
purely  algebraic  method  is  derived  in  the  frequency  domain  for  the  determina¬ 
tion  of  the  poles  and  transmission  zeros  of  a  matrix  transfer  function.  The 
matrix  Routh  algorithm  and  the  reverse  process  of  the  algorithm  (Shieh  and 
Gaudiano  1974,  Shieh  1975,  Shieh  et  al.  1975)  are  used  to  decompose  an  n0  x  ni 
rational  matrix  transfer  function  T(s)  into  Z),(s)_1AT,(s)  and  Nr(s)Dr{s)~l, 
where  the  polynomial  matrices  Dt(s)  and  AT,(s)  with  appropriate  size  arc  left 
co-prime  and  Nr(s)  and  Dr(s)  right  co-prime.  When  n0  =  n„  the  poles  (the 
transmission  zeros)  of  the  T(s)  are  determined  from  the  zeros  of  the  determinant 
D,(«)  or  Dr(s)  (AT,(s)  or  Nr(s)).  When  n0^n„  or  the  matrix  Routh  algorithm 
is  of  ill-conditioned  ease,  the  determinant  of  the  rectangular  polynomial 
matrices  N,(s)  and  Nr(s)  cannot  be  obtained.  An  n,  x  n„  matrix  transfer 
function  T+(s),  which  is  the  generalized  inverse  (Desoer  and  Schulman  1974) 
of  the  modified  T(s),  is  established  and  factored  into  Dl*(s)~'A' ,*(s)  when 
n0  ^  ni  or  Nr*(s)Dr*(s)~1  when  na  ^  The  transmission  zeros  of  the  T(s)  are 

determined  from  the  invariant  poles  of  the  T*  (,s)  or  from  the  zeros  of  the 
determinant  Dt*(s)  with  size  n,-  x  ni  or  Dr*(s)  having  size  n0  y  n0. 

Along  the  same  line,  recently,  several  approaches  have  been  proposed  by 
various  authors  (Kung  et  al.  1976,  Anderson  and  Jury  1976,  Emre  and 
Silverman  1976)  to  determine  the  relative  primeness  of  two  polynomial 
matrices.  The  generalized  resultant  matrix  (Barnett  1971)  and  the  generalized 
Bezoutian  and  Sylvester  matrices  are  used  in  their  works.  When  the  degree 
of  the  polynomial  matrices  that  might  have  a  high  degree  common  divisor  is 
high,  the  dimension  of  the  resultant  matrix  or  the  equivalent  test  matrix  is 
very  high.  As  a  result,  the  effectiveness  of  their  approaches  is  less. 

2.  The  matrix  Routh  algorithm  and  the  matrix  Routh  array 

In  a  single  variable  system  it  is  well  known  that  the  poles  and  zeros  of  a 
transfer  function  can  be  determined  from  the  respective  denominator  and 
numerator  polynomials  that  are  relatively  prime.  The  Routh  algorithm 
and  the  Routh  array  (Fryer  1959)  are  often  used  to  determine  the  common 
factor  of  the  two  polynomials  in  order  to  determine  the  pair  of  relatively 
prime  polynomials.  In  this  paper  we  extend  the  concept  to  a  multivariable 
system  that  is  described  by  a  matrix  transfer  function.  Let  us  define  that  R 
and  C  denote  the  field  of  real  numbers  and  complex  numbers,  respectively, 
and  i?[s]  and  R(s)  the  sets  of  all  polynomials  and  rational  functions  in  the  field 
of  complex  variables  having  real  coefficients.  We  also  define  that 
and  R(s)n »x"<  are  the  sets  of  all  n0xn,  matrices  with  elements  in  i?[s]  and 
R(s),  respectively. 

Consider  the  following  matrix  transfer  function  T(s)eR(s)n°xn<  which  is  a 
product  of  a  polynomial  matrix  A2(5)ei?ts]"'x"‘  and  the  inverse  of  another 
polynomial  matrix  A1(s)e^[s],x«,  where  q  =  min  (n0,  n,) : 

T(s)  =  At(s)  A^s)-1  =  [A21  +  A2ts  + ...  +  Ai%n8n~l] 

*[A,i  +  /li8«+...  +  AjtB+1«"]-1  (la) 
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T(s)  =  41(a)~142(s)  =  Mu  +  ^12s+ ...  +^,.n+,«B]-1 


x[^2,  +  ^22s  +  ...+^2iB«"-1]  (16) 

where 

h  n  + 1 

A2(s)=  X  ^2ii«i“1  and  ^4j(s)=  X 

i  =  i  i=i 


The  matrix  coefficients  in  the  ^t2(s)  and  At(s)  are  expressed  by  the  double 
subscript  notation  as  A2  IeJinoy"t  and  At  leR'IK'>  for  the  use  of  the  matrix 
Routh  algorithm.  If  the  T(s )  is  expressed  as  follows  : 


then 


and 


T(s)  = 


A0(s) 


O(s) 


A„(s)=  X  ^i(«)  =  I  V4'1  =  I 


n 


n 


0(S)=  X  <D^-1=  X 

i= 1  »=1 


(2) 


where  A0(s)e7?[s]  is  a  polynomial  and  is  an  identity  matrix.  By 

using  the  following  matrix  Routh  algorithm  and  the  reverse  process  of  the 
algorithm,  the  T(s)  can  be  factored  into  B,(s)-lNt(s)  and  A’r(s)Z)r(s)-1,  where 
Dt(s),  Nt(s),  Nr(s)  and  Dr(s)  are  polynomial  matrices  of  appropriate  size.  The 
matrix  Routh  algorithm  (Shieh  and  Gaudiano  1974)  and  its  reverse  process 
(Shieh  et  al.  1975)  of  the  algorithm  for  a  multivariable  system  (ni~nn)  are 
expressed  as  follows  : 


/^ll 

-^12 

-4 13  • 

•  -^l.n-^l.n+l 

Hi  ~  AiiA21  1  ( 

f 

\a 

21 

A  22 

-^23  ■■ 

■  A2n 

H2  —  A2lA3l  1  < 

f 

/^31  “  ^12  ™"  ^1*^22 

4  A  J 
**32  —  **  13 

^1*^23 

■^33  ■■ 

H3  =  A3lAn  l( 

f 

^41  ~  -^22  H 2-^32 

•^42  ~  *^23 

—  H  2^33 

(3  a) 

/^2n,l 

^Zn^  ^Zn,l^Zn+ltl  \ 

\  J 

"2n+l,l 

. 

The  Hj  in  eqn.  (3  a)  are  the  matrix  quotients.  The  block  elements  of  the  first 
and  second  rows  of  eqn.  (3  a)  are  the  matrix  coefficients  of  eqn.  (1  a).  The 
block  elements  of  the  subsequent  rows  are  evaluated  by  the  following  matrix 
Routh  algorithm  : 

Hi  =  AilAi+ for  i  =  l,2,...,2k  and  k^n  1 


rank  *4i+1>1  =n,- =  «0 


(3  6) 


Ai.j  —  Ai^j+i—Hi-  2-^i-i,>+i 


for  j=  1,2, 


*  =  3,  4, 


3  t  2 


When  the  two  matrix  polynomials  A^s)  and  At(s)  have  no  common  factor,  the 
matrix  Routh  array  will  terminate  normally  (i.e.  we  have  2 n  matrix  quotients). 
When  the  two  matrix  polynomials  have  a  common  factor  (the  common 
divisor),  the  matrix  Routh  array  in  eqn.  (3  a)  will  terminate  prematurely,  and 
the  last  non-vanishing  row  consists  of  the  matrix  coefficients  of  the  common 
factor  B(s)  in  the  original  matrix  polynomials  Ax(8)  and  A2(s).  If  we  have 
2k  matrix  quotients  Hjt  we  can  construct  a  pair  of  relatively  prime  matrix 
polynomials,  Nr(s)  and  Dr(s),  by  using  the  reverse  process  of  the  matrix  Routh 
algorithm  in  eqn.  (3  b) : 

Pm+i.i^I 

Pi,  i  =  HtPi+ul  for  t  =  2k,  2k  —  1 . 2,  1  (3  c) 

P i-2,j+i~ P \,i  +  B. i-iP for  i  =2fc+ 1,  2k . 3  ;  J  =  l,2 . fcj 

The  T(s)  in  eqn.  (1)  is 

T(s)  =  At(s)Ax(s)-' = Nr(8)B(8)[Dr(s)B(8)]-' = Nr(8)Dr(s)-' 

=  [Pt  i  +  P22S+...  +  Pt,iA>e~1][P\l  +  ■?,*«  +  ...+  (4) 

The  procedure  can  be  well  illustrated  by  the  following  numerical  example. 


Example  1 

Consider  that  the  common  divisor  B(s)  and  a  pair  of  relatively  prime 
matrix  polynomials  N,(s)  and  Dr(s)  of  the  following  matrix  transfer  function 
T(s)  are  required  : 

T(s)  =  At(8)A1(8)-l  =  NT(3)B(8)[Dr(8)BW]-^  =  Nr(8)Dr(s)~'  (5) 

where 

At(s)  =  An  +  Alzs  +  Al3s*  +  A^s3 

-(?  >0  -Xo  :5)'+(i  =0- 

=  -^21  ^22^  ^23^* 

=(.6,  =jm_\  -.vg 

n0=ni  =  2  and  n  =  3 


The  matrix  Routh  array  is 
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The  matrix  Routh  array  terminates  prematurely  because  the  only  one  block 
element  4gl  in  the  sixth  row  is  a  null  matrix  ;  therefore,  the  common  divisor 
B(s)  in  T(s)  is 

fi(s)  =  ^M  +  ^52s  =  (J  _\)  +  (J  ■{)*  (65) 

By  using  the  matrix  quotients  H1  ...  Hx  in  eqn.  (6  a)  and  applying  the  algorithm 
in  eqn.  (3  c)  we  have 

tff(s)  =  P21  +  /V  =  (_22  + 

and 

DM)  =  Pn  +  P„s  +  P13s2  =  (  j5  j  2)  +  (  - 1 

In  order  to  show  that  the  B(s)  in  eqn.  (6  6)  is  a  common  divisor  of  ^(s) 
and  in  eqn.  (5)  we  replace  AX  j  in  eqn.  (3  6)  by  Ptj  in  eqn.  (6  c)  and  apply 
the  algorithm  in  eqn.  (3  6)  to  eqn.  (6  c).  Thus  we  have  the  following  alternative 
matrix  Routh  array  that  has  the  same  matrix  quotients  H(  as  eqn.  (6  a)  : 


The  justification  for  the  B(s)  in  eqn.  (6  b),  which  is  the  common  divisor  of 
the  matrix  polynomials  in  eqn.  (5),  can  be  proved  by  the  following  induction 
method.  Since  the  matrix  Routh  algorithm  is  developed  from  the  repeated 
process  of  long  division  of  two  polynomial  matrices,  the  reverse  process  of  the 
algorithm  can  be  applied  to  eqns.  (6  a)  and  (6  d)  to  obtain  the  following 
identities  : 

./181  +  4S2s  =  P  61(.461  +  Asas)  =  P  gi-B(s)  =  B(s) 

An  +  =  Ht(As  j  +  45Js) 

=H,B(a) 

—  P  4i-®(®) 

+  -^32®  +  -^sa®2  ~D3(A  *i  +  At  2®)  +  +  -^M®) 

=  B3P 41 -B(®)  +  sP  slB(s) 

=  (^*31  +  P**®)P(®)  (®  e) 

Aqi  +  ^22®  +  ^23®2  =  D 2(^31  "b  -4s2®  +  -^33®*)  "b®(^41  +  ^42®) 

=  H3(P3i  +  P3is)B(s)  +  sP  4j  B(s) 

—  (P  21  +  -^>22®)-®(®) 

+Al2s  +  A 23s2  “f*  =  H i(A3 j  “b  -^22®  "b  -^23®2)  ~b  ®(-^31  "b  ^32®  *b 

=  H1(Pil  +  Pt3a)B(s)  +  s(P 31  + P  32s)B(s) 

—  (P11+P 12®  +  •f>is*2)-S(®) 

From  the  last  two  equations  in  eqn.  (6  e)  we  observe  that 

A2(s)^Nr(s)B(a)'j 

[  (6/) 

Al(a)^DT(S)B(8)) 

Therefore  B{s)  is  the  common  divisor  of  the  two  polynomial  matrices  .4,(s) 
and  A3(s). 

When  Wj#m0  and  rank  Ajtl¥=q  in  eqn.  (3  6),  the  matrix  Routh  algorithm  in 
eqn.  (3  b)  cannot  be  directly  applied.  The  matrix  Routh  algorithm  and  its 
reverse  process  of  the  algorithm  in  eqn.  (3)  are  modified  and  discussed  by  the 
following  case  studies. 

(1)  ^^(sMif®)-1  (7) 

where 


Case  1 


At(8)=  ^  i4s,<«<-1  and  At(8)  =  £  -^l.f®*-1 

i  “  1  >=1 


n0^n<,  T(8)eR(s)n«xn>,  ^42(s)e.ft[s]n«x,,<,  ^41(s)ei?[«]"<xn« 

T(s)  =  ^(sMi*®)-1  =  Nt(s)B(8)[Dr(s)B(8)}-'  =  Nr(s)Dr(s)-1 


where 


iVr(s)eif[«]»«x"<,  Nr(s)  =  £  ;  Z)r(s)e/f[s]»<x»- 

*“  1 

DM  =  *2  Pty~K  P1,k+i=In, ;  B(8)eR[8]'<x*<,  £(«)- "  %' Bf-* 
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The  B(s)  is  the  common  right  divisor  of  the  A2(s)  and  A,(s).  For  the  use  of  the 
matrix  Routh  algorithm,  the  matrix  coefficients  in  the  Nr(s)  and  are 

expressed  by  the  double-subscript  notation  as  P1 t  and  P2i  which  can  be 
obtained  by  the  algorithms  as  follows. 

The  matrix  Routh  algorithm  is 


Hi  =  AitiAi+1'l-1,  i=l,2 . 2 4  and  4«n 


rank  A(1  =n{ 


(8  6) 


j-1,2,  ....  t-3,  4,  ...J 

The  constant  matrices  Hi  with  appropriate  size  are  called  the  matrix  quotients. 
If  n0  >  »,-,  the  pseudo-inverse  of  =  =  is 

the  left  inverse  of  Ai+1 1. 

The  reverse  process  of  the  matrix  Routh  algorithm  is 


l  =  2k,2k- 1 . 2,1 


(8  c) 


P{— i,y+i  = %Pi— i  —  2k+l,  2k . 3,  j  —  1,  2,  k 


Case  2 

n0^no  T(s)eR(s)n»Xn‘,  i42(s)c^[s]*«xn»,  v41(s)eJ?[*]n<xn» 

T(s)  =  A^AM”'  =  D,(s)-'£(s)[N,(s)-'B(s)r'  =  2>t(«)-W,(s)  (9  a) 

where  t 

N,(s)eR[s]n'xni,  N,(s)  =  £  Qa.i8*'1’  D,(s)eB[s]n*xn* 

<= i 

!>,(*)=  ;  B{8)eB[s]n°xn°,  B(s)=  *£  Bjh* 

<=i  *=i 

For  the  use  of  the  matrix  Routh  algorithm,  the  matrix  coefficients  in  the  N ,(«) 
and  D,(s)  are  expressed  by  the  double-subscript  notation  as  Qt  i  and  Q2  i  which 
can  be  obtained  by  the  algorithms  as  follows. 

The  matrix  Routh  algorithm  in  eqn.  (8  6)  is  applied  to  determine  the  matrix 
quotients  Hi : 

Hi  =  AitlAMA~x,  <  =  1,2 . 2k  and  4<n  ] 


rank  ^441  =  n0 

=  —  J  =  ....  t  =  3,4,  ... 

The  new  reverse  algorithm  is 
Qtk+1.1  *^i»« 

«U  =  «»+i.iffn  1  =  24,2*-!,...,  2,1 


(9  6) 


(9  c) 


*‘Qi,f+Qi-i,i+iHi-v  *  =  24+  1,  24, ...,  3,  j  — 1,2,  ...,4 


(2) 

where 


(10) 


T(s)  =  A1(s)-'At(s) 


M  +  l 


42(a)=  £  and  4,(«)  =  £  A^-K 


t-l 


<=1 


Case  1 

n0 ^nit  T(8)eR(8)n*xni,  ^42(«)6i?[«]n»xns  41(a)e2?[*],l»xw® 

T(S)  =  ^1(S)-M2{5)  =  [B(«)Z)I(«)]-1fi(«)AT1(s)  =  2>l(®)-W,(*)  (Ha) 

where  k 

N,(s)eR[a]^i,N,(8)=  £  Qt  ^  ;  Z>,(«)e.R[«]B*XB» 

»-i 

dm-  *1  Qi.iS*-1,  Qi.k+i-ln, ;  2?(*)eRM-*B*,  £(«)= "  s'1^-1 

•=i  >=i 


The  matrix  coefficients  in  the  Dt(a)  and  N,(s)  are  expressed  by  the  double¬ 
subscript  notation  as  Q1§i  and  which  can  be  determined  by  the  following 
algorithms. 

The  new  matrix  Routh  algorithm  is 


Mi  =  A{+11~1  Aitl,  *  =  1,2 . 2Jfc  and  k^n 

rank  A{1  =  n0 


Ai,j  —  Ai_t'f+l~Ai-i'j+iMi_a,  j  — 1,2 .  t  — 3,  4,  ... 


(lift) 


The  constant  matrices  Mi  with  appropriate  size  are  called  the  matrix  quotients. 
If  n0<nt,  Ai+ttl~l  =  A{+t'tr[Ai+tlA(+lilT]-1  is  the  right  invesre  of  the  Ai+tl. 
The  reverse  algorithm  in  (9  c)  is  applied  to  determine  the  Qlfi  and  Q2  i. 


^2*+l,l—^»o 

Qi,i  =  Qi+i,i^it  1  =  21,  2i-I, ...,  2, 1 

i  —  2k+  1,  2k, ...,  3,  j  =  1,2 . k J 


(11c) 


Case  2 

«<,>»<,  2,(a)e-R(a)B»XB<,  ^4,(s)eR[a]B«XB',  ^41(a)eR[a]B«XB» 

T(s)  -  AM-1  AM  =  [R(a)iVr(a)-‘]-‘R(a)2>r(a)->  =Nr(s)Dr(s)^  (12  a) 

where  t 

iVr(a)£R[«]B»XB<,  Nr(s)  =  j  P^a*-1;  Z?r(a)eR[a]"<XB' 

<-i 

aw=  s'  v1.  Pi.k+i=in,  ->  mem***”,,  B(a)= 

i-l  »«1 

The  matrix  coefficients  in  the  Z?r(«)  and  Afr(«)  are  expressed  by  the  double¬ 
subscript  notation  as  PlA  and  P"  which  can  be  obtained  by  the  algorithms 
as  follows. 
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The  matrix  Routh  algorithm  in  eqn.  (116)  is  applied  to  determine  the 
matrix  quotients  : 

Mi  =  Ai+li~1Ail,  *  =  1,2 . 2k  and  k^n 

rank^li>1=»1  (12  6) 

■Ai,j  =  Ai_tj+1  —  Ai^lii+lM{_t,  j  =  l,2 .  *  =  3,4,... 

The  reverse  algorithm  in  eqn.  (8  c)  is  applied  to  determine  the  Pl  f  and  Pt  i : 

Pik+l.l  =  Im 

l  =  2k,2k-l . 2,1  (12  c) 

P i-jj+i  =  P jj  +  i-i,/+i,  *  =  2t+l,  2i,  3,  .7  =  1,2 . fc 

By  using  Gilbert’s  theorem  it  has  been  shown  (Shieh  and  Gaudiano  1975)  that 
the  dynamic  state  equations,  which  are  constructed  by  using  2k  matrix 
quotients  H{  or  Mt  that  are  obtained  from  the  matrix  Routh  algorithms, 
are  minimal  realizations  of  the  T(s).  The  minimal  dimension  of  the  system 
matrix  is  kq  x  kq,  where  q  =  min  (nf,  n0)  and  kq  =  rank  T(s)^r0.  The  rank  T(s) 
can  be  determined  from  the  Hankel  matrix  (Ho  and  Kalman  1966).  By 
using  the  same  2k  matrix  quotients  H{  or  Mt  and  performing  the  reverse 
process  of  the  matrix  Routh  algorithm,  the  monic  polynomial  matrices  Dr(s) 
and  /),(«)  are  obtained  and  shown  in  eqns.  (8),  (9),  (11)  and  (12).  The  highest 
power  of  8  in  the  Dr(s)  and  D,(s)  is  k  and  the  matrix  coefficients  of  sk  (i.e. 
Plwk+ 1  and  Qiik+i)  are  identity  matrices  having  size  qxq.  Therefore 

kq  + 1 

det  Dr(s)  =  det  D,(s)  =  £  d^-'ldfs)  (13) 

i-i 


The  highest  power  of  s  in  the  monic  polynomial  d(s)  is  kq,  which  is  the  rank  of 
the  T(s).  As  a  result,  the  d(s)  in  eqn.  (13)  is  the  characteristic  polynomial  of 
the  T(s)  and  the  polynomial  matrices  Dr(s)  and  Nr(s)  are  right  co-prime  and 
the  D,(s )  and  N,(s)  are  left  co-prime. 

From  the  above  discussion  we  also  note  that  the  necessary  condition  for  the 
existence  of  the  matrix  Routh  algorithm  is  that  the  ratio  (denoted  as  k)  of  the 
rank  T(s)  and  the  minimal  dimension  of  the  T(s),  rjq=k,  is  an  integer.  If 
the  ratio  r0/q  is  not  an  integer  or  it  is  an  integer  but  the  condition  (rank  A(l  =  q) 
in  the  matrix  Routh  algorithm  in  eqns.  (8)— (12)  is  violated  due  to  the  ill- 
conditioned  matrix  Ai  v  then  the  original  T(s)  should  be  modified.  T(s)  is 
modified  by  adding  another  transfer-function  matrix  T0(s)  =  ljg(s)K  whose 
rank  is  of  (kq  —  r0),  where  k  is  the  nearest  integer  and  the  scalar  polynomial 
g(s)  is  not  a  factor  of  the  A0(s)  in  eqn.  (2).  The  A  is  a  constant  matrix  with 
appropriate  dimension.  The  modified  system  T1(s)eB(s)n«xni  is 

THs)  =  T(s)  +  ±K  (14) 

where  rank  [(l/gr(s))ZL]  =kq-r0  and  rank  T1(s)=kq.  It  is  noted  that  the 
addition  of  (1  jg(s))K  to  the  T(s)  does  not  affect  the  locations  of  the  poles  of  the 
T(s). 
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3.  Determination  of  poles  and  transmission  zeros 

By  using  the  matrix  Routh  algorithm  the  T(a)  is  factored  into 
and  N r(s)Df(s)-1,  where  D,(a)  and  AT((s)  are  left  co-prime  and  Nr(a)  and  Dr(s) 
are  right  co-prime.  When  =  n0  =  q,  Desoer  and  Schulman  (1974)  have  shown 

that  the  transmission  zeros  of  the  T(s)  are  the  zeros  of  the  scalar  polynomial 

n^'  °r  n(s)  =  det  Nt(s)  =  det  Nr(a)  =  0  (15) 

where  Nt(a),  Nr(a),  D,(a)  and  Dr(a)  are  polynomial  matrices ;  N,(s)  and  Nr(s) 
are  q  x  q  ;  D,(«)  and  Dr(a)  are  qxq.  The  poles  of  the  T(a)  are  the  zeros  of  the 
following  characteristic  equation  : 

A  (a)  =  det  D,(a)  =  det  Dr(a)  =  0  (16) 

When  rjq^k  (an  integer),  the  matrix  Routh  algorithm  cannot  be  applied. 
The  procedure  shown  in  eqn.  (14)  can  be  applied  to  determine  a  pair  of  relatively 
prime  polynomial  matrices  Dtl(a)  and  Nfta)  or  Nrl(a)  and  Drl(a)  as  follows  : 

T\s)  =  =  Nrl(8)Dr'(a)-'  (17  o) 

The  poles  of  the  T^a)  can  be  determined  from  the  following  equations  : 

det  /),*(*)  =  det  Dr1{s)  =  {g(a)}lal~T»P(a)  =0  (17  6) 

where  the  g(s)  is  the  polynomial  used  in  eqn.  (14).  The  poles  of  the  T(s)  are 
the  zeros  of  P(a)~0. 

When  rjq  =  k  (an  integer)  and  n0/«<(  the  Nt(a)  and  Nr(a)  obtained  from  the 
matrix  Routh  algorithm  are  not  square  polynomial  matrices  of  size  n0  x  n0  and 
ni  x  rif.  Therefore  the  transmission  zeros  cannot  be  directly  determined  from 
eqn.  (15).  The  transmission  zeros  of  the  T(a)  can  be  determined  from  the 
invariant  zeros  of  the  determinants  of  all  q  x  q  minors  of  the  N,(a)  or  Nr(a)  in 
eqn.  ( 1 5)  where  q  =  min  (n0,  n4).  However,  when  the  nt  is  much  larger  than  the 
w0  and  vice  versa,  there  exist  many  qxq  minors  which  are  expressed  by  poly¬ 
nomial  matrices. 

It  is  a  cumbersome  task  to  find  the  determinants  of  these  qxq  minors  and 
to  determine  the  roots  of  many  polynomials.  This  difficulty  can  be  overcome 
by  applying  Desoer  and  Schulman ’s  (1974)  theorem.  The  transmission  zeros 
of  the  T{a)  are  obtained  from  the  invariant  poles  of  two  generalized  inverses 
of  the  modified  T[a).  In  this  paper  we  present  a  procedure  to  obtain  the 
generalized  inverses  of  the  modified  T(a)  so  that  the  transmission  zeros  of  the 
T(a)  can  be  determined.  The  steps  are  shown  as  follows. 

Step  1.  Modify  the  T(a)  and  formulate  the  generalized  inverses  Tt*(a), 
i  =  l,  2  of  the  modified  T(a),  or 

r1*(a)=K(s)r(«)]-»  =  41(«)tm1(s)41(a)]-1,  »„»»,  (18  a) 

=  tw<(#)r(s)]-,  =  tTO<(«)^1(s)]-M1(«),  n0<«<  (18  6) 

where  Ti*(a)eR(a)n*^n*,  -d,(s)e.R[a]«*'',  /fa(«)eU[s]B*x"<,  m{(«)eJR[s].  The  monic 
polynomials  mi(«)€i{[«],  *  =  1,2,  should  not  be  the  factors  of  the  \(a)  in  eqn.  (2) 
and  of  the  g(a)  in  eqn.  (14),  but  are  chosen  in  such  a  way  that  the  power  of  the 


i 

L 
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s  in  the  polynomial  matrices  mi(a)Ai[s)  in  eqn.  (18)  is  1°  higher  than  that  of 
the  At(«).  This  modification  does  not  affect  the  transmission  zeros  of  the 
T(8)  but  it  makes  the  matrix  Routh  algorithm  applicable. 

Step  2.  Modify  the  Ti*(s)  when  rjq^k  (an  integer). 

The  process  of  this  modification  is  shown  in  eqn.  (14).  The  additions  of 
(1  lg(s))K,  i  =  1,  2  to  the  T+ls),  i  =  1,  2,  do  not  affect  the  locations  of  the  poles 
of  the  T4*(«)  or  the  transmission  zeros  of  the  T(s).  After  performing  some 
matrix  operations  we  have  the  generalized  inverses  of  the  modified  T(s) 
denoted  as  Ti+(s),  *  =  1,2: 

r<+(«)  =  [?j(*)^i(®)  +  wii(«)^^»(*)][?.(s)wt<(*)^»(*)]-1.  n0>n(  (19  a) 

Step  3.  Determine  two  pairs  of  relatively  prime  polynomial  matrices. 

The  algorithms  in  eqns.  (9)  and  (12)  can  be  applied  to  obtain  two  pairs  of 
left  co-prime  polynomial  matrices  denoted  as  Du*(s)  and  Nu*(s),  *  =  1,2  or 


right  co-prime  polynomial  matrices  Nri*(s)  and  !>**(«),  *'  =  1,2: 

Ts+(«)  =  Dli*(8)~1N  h*(8),  n02nt  (20  a) 

=NH*(s)Dr(*(8)-1,  n0Hn{  (20  6) 

Step  4.  Select  the  required  transmission  zeros. 

The  poles  of  the  Tt+(s),  *  =  1,2,  are 

det  Dh*(8)  =  (m<(«))«n(s)i:i(s)  =  0,  »0  >  n{  (21  a) 

det  Dh*(«)  =  (w1(s)}9n(s)fcj(s)  =  0,  n0  ^  nt  (216) 

The  required  transmission  zeros  of  the  T(s)  are  the  invariant  poles  of  the 
Tt+(s),  *  =  1,2.  They  are  the  zeros  of  polynomial  n(s)  in  eqn.  (21),  or 

n(s)  =0  (21c) 


When  the  T(s)  is  not  a  strictly  proper  rational  matrix  transfer  function,  eqn. 
(19)  can  also  be  used  to  determine  the  poles  and  transmission  zeros  of  the  T(a). 


Example  2 

Consider  that  the  poles  and  transmission  zeros  of  the  following  matrix 
transfer  function  T(s)  are  required  to  be  determined  : 

r(s)  =  As(s)A1(s)-‘  (22) 

where 

IV  +  6«*+ll*+12  «*  +  9s  +  20  1 


A,(«)  = 


«*  +  4«-6  s*  +  3«*-7a  +  15 


s*  +  6s*+13s+10  3s*  +  6«  +  26 
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and 


s4  +  5s3+lls*+13s  +  6 
s3  +  2s*  —  s  —  2 


s3  +  8s*  + 17s  +  10 
8*  +  6s  —  s2  +  5s  +  6 


n0  =  3  and  nt  =  2 


It  is  difficult  to  apply  most  time-domain  approaches  to  this  problem.  By 
using  the  proposed  algorithms  in  eqn.  (8),  the  T(a)  can  be  factored  into  a  pair 
of  relatively  prime  polynomial  matrices  A?r(s)  and  Dr(s) : 


T(s)  =  =  Nt(8)B(a)[Dr(8)B(8)}-'  =  A»Z>,(s)->  (23) 

where 


"s  +  4 

0 

"s*  +  3s  +  2 

0 

A»  = 

0 

s  +  4 

s  +  5 

2 

and  Dr(a)  = 

0 

s*  +  3s  +  2 

Following  eqn.  (16)  we  have  the  required  poles  of  the  T(s) : 

A (s)  =  det  Dr{a)  =  (a  + 1  )s(s  +  2)*  =  0  (24  o) 

or 

Sj=s8=-1  and  s3=e4=-2  (24  6) 


It  is  interesting  to  note  that  the  common  factor  B(a)  in  eqn.  (23)  is 

fs4  +  2s  +  3  s  +  5  T 


B(s)  = 


8—1 


a3  —  2s  —  2 


(25) 


Since  n0  >  ni  and  the  determination  of  the  transmission  zeros  of  the  T(a)  are 
required,  we  construct  the  generalized  inverses  T{*(a)  in  eqn.  (18)  from  the 
modified  T(s)  in  eqn.  (23),  and  apply  the  proposed  algorithm  in  eqn.  (9)  to 
decompose  the  Tj*(s)  into  two  pairs  of  left  co-prime  polynomial  matrices 
Du*(a)  and  NH*{8),  i  =  1,  2  in  eqn.  (20  a)  as  follows  : 


Tx*(a)  -  ^(sWaM.fs)]-1  =  =  Dn*{s)~'N  n*(s)  (26  a) 


where  mt(a)  =  s3  +  s  + 1  is  not  a  factor  of  the  A  (a)  in  eqn.  (24  o) : 

5s2  +  5s  +  4  l-I965s2  +  2-1649s  +  0 

0  s3  +  5-88073s*  +  5-88073s 

329s +  0-981  0-01776+0 
1407s -0-574  s*  +  3-04s  + 


[s3  +  i 
D,S(s)=\ 

[0-4275s!  + 1-3 
00797s1  +  0-1 


•81525  ] 

+  4-88073 J 
■057  0-573s*+ 1-67  Is +1-02  "1 

1-72  -  0-08s* -  0- 1 43s  +  0-574 J 


Tt*(a)  =  ^(stfw^sM ,(s)]-»  =  Dr[8)[mt[a)N r(s)]_1  =  Dlt*(a)-'Nn*(a)  (26  6) 
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where  w2(s)  =s2  —  10  is  not  a  factor  of  the  A (s)  in  eqn.  (24  a)  : 


£>«*(«) 


-s3  +  4s2  _  10*  _40  -  3-854672«2  +  38-54672 

0  s3  +  0-4438694s2  —  10s  — 4-438694 


Nl2*(s) 


"2-24s2  +  4-45s  +  1-41  -1-38S-1-78  -  1 -24s2- 1 -45a +  0-59 

1  -02s2 -0-36s- 0-04  s2  +  0-49a  +  0-16  -  l-02s2 +  0-36* +  0-038 


Following  eqn.  (21  a)  yields 

det  Dn*(s)  =  (a2  +  a  +  1  )2(a  +  4)(a  +  4-88073)  =  w,2(s)«(a)/fc,(a)  =  0  (27  a) 

det  Z>,2*(a)  =  (a2  -  10)2(a  +  4)(a  +  0-44386943)  =  m22(a)n(a)A;2(a)  =  0  (27  b ) 

The  common  divisor  n(a)  of  the  det  Z>„*(s)  and  det  D,2*(s)  is  the  common 
factor  «(a)e/?[a]  in  eqns.  (27  a)  and  (27  b ).  The  transmission  zeros  of  the  T(s) 
which  are  the  invariant  poles  of  the  31i*(a)  in  eqn.  (26)  are  the  zeros  of  the 
polynomial  n(a),  or 

n(a)  =  a  +  4  =  0  (28) 

The  transmission  zero  is  s  =  —  4. 

The  computation  involves  only  arithmetic  operations  of  small-size  matrices. 
Therefore,  it  is  believed  that  the  proposed  method  is  computationally  superior 
to  most  time-domain  approaches  if  the  system  is  given  in  the  frequency  domain. 


4.  Conclusion 

A  purely  algebraic  method  has  been  presented  for  factorizing  a  rational 
matrix  transfer  function  into  a  pair  of  relatively  prime  polynomial  matrices 
and  for  determining  the  poles  and  transmission  zeros  of  a  multivariable  system. 
Also,  the  common  divisor  of  two  matrix  polynomials  can  be  determined  from 
the  matrix  Routh  algorithm  and  the  matrix  Routh  array.  When  a  matrix 
transfer  function  that  might  have  a  high  degree  common  divisor  is  given,  the 
method  proposed  in  this  paper  is  computationally  superior  to  most  time- 
domain  methods  because  the  proposed  algorithm  only  deals  with  arithmetic 
operations  of  small-size  matrices.  The  matrix  Routh  algorithm  has  been 
extended  for  general  cases  (w,-  #  n0  and  n(  =  n0). 
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Some  sufficient  and  some  necessary  conditions  for  the  stability  of  a  class  of  multi- 
variable  systems  represented  by  matrix  polynomials  are  derived.  .4  new  linear  block 
transformation  is  also  established  for  transforming  an  observable  block  companion 
form  to  the  block  Schwarz  form. 


I  Introduction 

The  accurate  description  of  most  practical  systems,  for  ex¬ 
ample  both  a  small  semiactive  terminal  homing  missile  system 
(IJ  and  an  aircraft  system  [21,  result  in  high  order  coupled  mult  i- 
variable  differential  equations.  Linear  representations  of  these 
systems  are  by  a  set  of  coupled  high-order  differential  equations 
or  a  matrix  differential  equation.  A  primary  concern  in  the  design 
of  these  multivariable  systems  is  the  stability  problem.  One  con¬ 
ventional  approach  is  to  formulate  the  system  into  a  high  dimen¬ 
sional  state  equation,  then  to  determine  the  stability  by  either 
directly  evaluating  the  roots  of  the  scalar  characteristic  poly¬ 
nomial,  indirectly  applying  the  Uouth  criterion  [3],  or  applica¬ 
tion  of  Jury’s  inner  theory  [4]  on  the  characteristic  polynomial. 
However,  the  determination  of  a  characteristic  polynomial  for 
a  large  dimensional  system  is  tedious.  Moreover,  if  a  system  is 
modeled  as  a  matrix  differential  equation,  it  is  more  natural  to 
determine  the  stability  directly  from  the  matrix  polynomial  than 
indirectly  from  a  scalar  polynomial,  riome  approaches  have  been 
proposed  to  determine  the  stability  of  a  multivariable  system 
directly  from  the  matrix  polynomial.  Papaconstantinou  [5| 
suggested  a  scheme  for  testing  stability  of  polynomial  matrices. 
In  his  work,  a  recursive  algorithm  was  developed  to  compare  the 
normalized  largest  eigenvalues  with  unity.  However,  the  method 
requires  the  calculation  of  the  eigenvalues  of  largest  moduli  for 
indirectly  determining  the  stability  of  polynomial  matrices.  Re¬ 
cently,  Shieh  and  Sacheti  [6]  partially  extended  the  scalar  lloulh 
criterion  [3]  to  the  matrix  case.  In  this  work,  it  is  shown  that,  if 
a  matrix  polynomial  B(s)  =  Is"  +  Bns"~l  +  . . .  -j-  li,  is  given,  a 
matrix  Routh  array  can  be  constructed  by  using  the  following 
matrix  Routh  algorithm: 

Ci.i  **  B.+ t-n  j  ”  1,  2,  3,  . . .,  I 
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J+» 

n  even 

where  l  = 

In  -f  1 

[  2 

n  odd 

L'i.y 

—  iin+t-ij  j  ~  3,  .  . 

1 

Cn 

=  I 

Ci.i 

i  «  1,2,.. 

•  .  i  =  3. 

4,  ... 

=  CuCCW,.!)-!  t  -  1,  2, 

...,  n 

det  (C,+ ],i)  5*  0  (1) 

A  sufficient  condition  for  stability  of  the  det  [B(s)J  is  that  all 
tire  “matrix  quotients"  //,  be  real,  symmetric,  positive  definite 
matrices.  Note  that  this  sufficient  condition  deals  only  with  H, 
and  not  (the  block  elements  in  the  first  column  of  the  matrix 
Routh  array).  Liapunov  theory  with  the  state  equation  in  the 
controllable  block  companion  (controllable  phase-variable)  form 
was  used  to  derive  their  result. 

In  this  paper,  we  develop  two  approaches  for  determining  the 
stability  of  a  class  of  multivariable  systems.  One  approach  uses 
the  “matrix  quotients”  .1/,  t hat  arc  developed  from  an  alternate 
matrix  Routh  algorithm  and  a  state  equation  in  the  observable 
block  companion  form  (7|.  The  other  approach  uses  the  block  ele¬ 
ments  in  the  first  column  of  the  matrix  Routh  array.  Two  suf¬ 
ficient  conditions  and  three  necessary  conditions  are  derived 
for  the  stability  of  matrix  polynomials,  thereby  partially  ex¬ 
tending  the  scalar  Routh  criterion  to  the  matrix  Routh  criterion. 

II  Sufficient  Conditions 

The  objective  of  this  paper  is  to  establish  the  criteria  for  the 
stability  of  the  following  matrix  differential  equations. 

•+i 

£  -  [01,  B,»  -  I  (2a) 

i-1 

Transactions  of  the  ASME 


and 


f),_1x(0)  =  [«i_ 


i  =  1,  2,  3,  ....  n 


(26) 


where  x(f)  is  the  m-dimensional  state  vector.  /!,,  1,  and  |0)  are 
in  X  m  real  constant  niatrix,  identity  matrix  and  null  matrix, 
respectively.  For  the  scalar  case,  it  is  well  known  that  a  system 
is  asymptotically  stable  if  and  only  if  the  Routh  array  elements 
in  the  first  column  are  all  positive.  Shieh  and  Sacheti  [<>]  partially 
extended  the  Routh  criteria  (31  to  the  matrix  case  and  derived 
a  sufficient  condition  for  the  stability  of  a  multivariable  system 
in  equation  (2)  from  the  controllable  block  companion  form.  In 
this  paper  we  derive  some  sufficient  and  some  necessary  condi¬ 
tions  for  the  system  in  equation  (2)  from  the  observable  block 
companion  form. 

Let  us  rewrite  the  system  in  equation  (2)  into  the  following 
observable  block  companion  form: 


algroithm  and  alternate  matrix  Routh  array  which  are  different 
from  those  in  equation  (1). 

Let  us  define  I  —  n/2  +  1  if  «  is  an  even  number,  otherwise 
I  =  (a  +  1  )/2,  and  l),.,  as  follows: 

«!.,  =  j  =  1,  2,  3,  ....  I 

D,.i  =  )  =  1,  2,  3,  ...,; 

Du  =  1  (5a) 

The  alternate  matrix  Routh  array  and  the  matrix  Routh  algo¬ 
rithm  are: 


-Du  “ 


N  -  [BIN 


(3a) 


Mi  —  < 


(2(0)1  - 

M 

(36) 

Dn  -  B. 

where 

M,  = 

Di.-tthi  < 

jDh  ^  Dn  —  OjtAf  x 

0  0  0 

■  0  ^ 

I  0  0 

■  0  —B\ 

Mi  - 

D*rlDn  < 

IB]  = 

0  10 

■  0  -B, 

Mi  - 

Da-'D„  < 

0.40,-  DnMi 

0  0  0 

■  I  -B„_ 

Dh  ^  Da  —  DoMi 

The  dimensions  of  the  matrix  [B] 

the  block  elements  B„  and 

state  vector  [x]  are  ( nm )  X  (nm), 

m  X  m,  and  (nm)  X  1,  re- 

Du  i 

spectively. 

Equation  (3)  can  be 

transformed  into  the  block 

Schwarz  form  by  using  the  following  linear  transformation : 

Af.  - 

Dn+i,rl^n.X 

< 

N  =  [*,][»] 

(4ay 

A»+i,i 

and 

where 

[»]  = 

[tf.l-MBHtfdl,,]  -  [Ally] 

(46) 

where 

I 

0  1 

0 

0 

i 

DuDa~l 

0  Di)Dn~l 

0  DuD„- 1 

0 

0 

0 

DuDn~l 

0 

DnD,r'  0  | 

[if,)  = 

0 

0 

I 

0  DnD4rl 

0  DuDu'1 

0 

0 

0 

:  i 

0 

Z>wD.r»  0 

0 

• 

0 

0 

1  0  ! 

I 

1 

0  DnDiC1 

0 

• 

0 

0 

0  ! 

0 

t 

1 

1 

I  0 

0 

0 

0 

:  o  ; 

0 

l 

1 

0  !  I 

On  =  Bh~x 


Dn  -  B. 


hu  —  B*-t. . . 


Dn  m  B*~i  . . 


(56) 


(4c) 


and 


0  -di  0-0  0 

I  0  -A,  ■  0  0 


0  0  0  0  -/I*., 

0  0  0  •  I  -A, 


The  dimension  of  each  block  element  in  (A]  and  \K,)  is  m  X  m. 
The  block  element*  having  dimension  m  X  m,  in  c(|iintion 
(4c)  can  be  obtained  from  the  following  alternate  matrix  Routh 


D<j  —  D,-,.y+1  —  j  •  1,  2,  . . .,  t  —  3,  4,  . . . 

Mi  —  Di+i,i_1D(,i  »  =  1,  2,  . . .,  n 

det  [Divi.il  *  0  (5c) 

The  construction  of  the  matrix  Routh  array  in  equation  (56)  is 
as  follows.  Arrange  the  matrix  coefficients  of  the  given  matrix 
polynomial  in  equation  (2a)  in  the  first  two  rows  of  the  array 
shown  in  equation  (56).  A  new  matrix  Mi  is  obtained  by  the 
matrix  multiplication  Dir'On  where  Du  and  Du  are  the  block 
elements  in  the  first  column  of  the  array.  The  block  elements  in 
the  third  row  are  generated  from  the  Mi  and  the  block  elements 
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in  tKe  first  two  rows  as  follows:  First,  each  block  element  in  the 
second  row  is  postmultiplied  by  .l/i,  then,  subtract  each  resulting 
matrix  from  each  block  element  in  the  first  row,  finally,  shift 
each  block  element  so  obtained  one  column  left  and  drop  the 
zero-first  block  element  to  form  the  third  row.  The  second  and 
the  obtained  third  row  are  then  used  as  starting  rows  to  generate 
the  new  matrix  .W,  and  the  oloek  elements  in  the  fourth  row.  Re¬ 
peating  the  processes  to  the  n  1  row  yields  the  complete  matrix 
l’outh  array.  When  any  matrices  other  than  On  or  0„ti,i 

in  equation  ( -V )  are  singular,  another  set  of  11 .11  can  be  chosen 
from  the  new  matrix  polynomial  that  is  the  product  of  the 
original  matrix  polynomial  and  att  asymptotically  stable  matrix 
polynomial.  Thus  a  new  matrix  liouth  array  can  be  obtained 
and  the  stability  of  the  original  matrix  polynomial  is  preserved 
because  the  stability  of  the  original  matrix  polynomial  is  in¬ 
variant  under  this  transformation.  Shieh  and  Sacheti  [6]  have 
shown  that  if  //,  =  C,,iC\+m~1  for  »'  =  1,  2,  ....  n  in  equation 
(1)  are  positive  definite,  then  the  system  in  equation  (2)  is  stable. 

Here,  we  show  similar  results  when  replacing  Hi  by  Mi.  Note 
that  a  positive  definite  matrix  means  a  matrix  is  real,  symmetric 
and  positive  definite. 

Ttiaoram  1.  If  { M i  =  1,  2,  . . .,  n  in  equation  (5)  are  posi¬ 
tive  definite,  then  the  system  in  equation  (2)  is  stable. 

Proof.  Performing  the  following  new  transformation 

(y)  -  l*.]|z]  («) 

on  equation  (4)  yields 

M  =  (fc,]-i[d][/r,)[*i 

=  [Fl!*]  (7a) 


where 


i  0  0  0  ' 

0  D„-,,i  •  0  0 

[*.)  - .  (76) 

0  0  Du  0 

0  0  0  2),, 


and 


where 

~M.  0  •  0  0  " 

0  Af._,  •  0  0 

P- .  (86) 

0  0  ■  M,  0 

0  0  •  0  M, 

and  T  in  equation  (8a)  designates  transpose. 

Since  |  Af  ,|  are  positive  definite  which  implies  that  P  is  posi¬ 
tive  definite,  V  is  positive  definite.  The  derivative  of  V  is 

V  =  \x]T[PF  +  F'7>][z] 

=  -  MHOllr]  =  -  Mr[ft«r][*]  (9a) 

where 

'o  o  •  o  o  ~i  r°" 

0  0-0  0  0 

IQ1  . . ,  \R)  ■=  •  (96) 

0  0-0  0  0 

0  0-0  2/J  V2/_ 

rank  [Q]  =  rank  [/?]  =  m.  From  equations  (8)  and  (9)  we  can 
see  that  V  is  a  Liapunov  function.  Hence,  we  conclude  that  the 
system  in  equation  (2)  is  stable. 

From  the  result  obtained  in  Theorem  1,  we  establish  another 
sufficient  condition  for  the  stability  of  the  system  in  equation 
(2)  by  using  the  block  elements  /),.,  in  the  matrix  Routb  array 
in  equation  (5)  instead  of  the  Mi  in  equation  (5). 

Theorem  2.  If  {Z>i.,|  i  =  2,  4, 6,  . . .,  are  positive  definite,  the 
eigenvalues  of  { D.i)  i  =  1,  3,  5,  . . .,  are  positive  and  real,  and 
) D(.iD>+i,i)  t  =  1,  3,  5,  . . {2),+i.)Di.i|  i  =  2,  4,  6,  . . are 
Hermitian,  the  system  [equation  (2)]  is  stable. 

In  order  to  prove  Theorem  2,  we  need  the  following  lemma 
which  is  due  to  KyFan  [9]  [p.  137). 

Lemma  X.  Let  Ki  be  positive  definite  and  Kt  such  that  KtKi 
is  Hermitian.  Then  K,K»  is  positive  definite  if  and  only  if  the 
eigenvalues  of  A,  are  positive  and  real.  In  the  following  lemma, 
we  switch  the  conditions  on  A,  and  Kt  yielding  the  same  result. 


It  is  noticed  that,  if  each  block  element  in  the  matri.:  [F[  in 
equation  (7c)  were  a  scalar,  then  the  matrix  [F]  would  be  a 
matrix  of  the  Schwarz  form  [S|,  Since  the  elements  are  blocks, 
the  matrix  |F]  in  equation  (7c)  is  a  block  Schwarz  form  matrix. 

The  linear  transformation  matrix  [A)  between  x  coordinates 
and  z  coordinates  is 

|r|  -  [A|[zJ  -  ; ATil[ (7rf) 
Now,  consider  the  following  quadratic  equation : 

V  -  [z]r[F][z|  (8o) 


Lemma  2.  Let  Kt  be  positive  definite  and  FT,  such  that  it, A', 
is  Hermitian.  The  A,K ,  is  positive  definite  if  and  only  if  the 
eigenvalues  of  Ki  are  positive  and  real. 

Pi— t.  Since  Kt  Is  positive  definite  which  implies  KtT  is  posi¬ 
tive  definite,  where  T  designates  transpose,  it  is  seen  from  lemma 
1  that  A'irA',r  is  positive  definite  if  and  only  if  the  eigenvalues 
of  A'ir  are  positive  and  real.  But  A,rA'ir  (A',A',)r;  i.e., 
KiKt  is  positive  definite  if  and  only  if  the  eigenvalues  of  Ki  are 
positive  and  real. 
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Lemma  3.  If  A'j  is  positive  definite  and  A'iA'i  is  symmetric, 
then  A'r'A',  i-s  symmetric. 

Proof.  Since  (KtK,)T  =  K.TKJ  =  K.KJ  =  A’lAj  which 
implies  KiT  =  A'j  1  A', A..  Hence  (At  *A'j)r  --  K»T(Ki  *)r 
=  A'jA'r1  A’l1  A'j  =  A't_lA'j;  i.e.,  A'r'A'j  is  symmetric. 

Proof  of  Thaoram  I.  By  lemma  3,  we  know  that  D;*.i.rlD,.i 
is  symmetric  for  i  =  1,  2,  ....  By  lemma  2  or  3,  we  know  that 
.If,  =  fliji.r'Oi.i  is  positive  definite.  Hence,  the  system  in 
equation  (2)  is  stable  following  the  results  of  Theorem  1 . 

In  order  to  show  an  application  of  Theorem  1  and  Theorem  2, 
let  us  consider  the  following  matrix  characteristic  equation: 

A&  +  Us  +  C  =  (I  (!<*»> 


-  -  C  "]■  *  ■  C‘  J]  "d  c  ■ 


If  we  arrange  the  matrices  .4,  Ht  and  C  in  equation  (U)fl)  by 
following  the  matrix  Itouth  algorithm  of  equation  (1),  we  obtain 


-i[-r 


c„  -  c 


r 

.0  2. 


In  Bellman  [9],  Ip.  f»7,  p.  1011  it  is  shown  that,  if  At,  Bt,  and  C\ 
arc  positive  definite,  then  the  roots  of  det  [AiS1  +  Bis  +  Ci)  =  0 
have  negative  real  parts.  But  in  this  example,  no  conclusion  can 
be  made  from  Bellman’s  results.  However,  we  know  that 

„  „  ,  „  P 25.77  13.71 

Ou  •  =*  At  •  Bt  —  I  I 

|_13.7  7.3J 

and 

rs  1.631 

D«  •  Du  -  C,  ■  B,  =  (116) 

1.1.63  0.9  J 

which  are  symmetric,  Bt  is  positive  definite,  and  the  eigenvalues 
of  .-1i  and  f'i  are  posit  ivc  and  real.  Therefore,  from  Theorem  2 
we  conclude  that  the  system  in  equation  ( 1  la)  is  stable.  Although 
only  second  order  matrix  polynomials  with  2X2  matrix  coef¬ 
ficients  are  illustrated  in  the  examples,  the  theory  is  valid  for 
high  order  matrix  polynomials. 


p  °1 

O 

il 

O 

II 

r  n 

1  3 

Lo  U 

_0  2. 

r* r 

Li  3. 

In  this  case,  no  conclusion  can  be  drawn  from  the  sufficient  con¬ 
dition  established  by  Shieh  and  Sacheti  [61.  However,  if  we  ar¬ 
ranged  the  matrices  A,  B,  and  C  according  to  equation  (5),  we 
have 


[o  J  D»~c~  1  3 

LO  2. 


A/,-- 


D«  =  D  =  [J  ‘] 

i] 


From  Theorem  1,  we  see  that  the  system  is  stable. 

This  example  shows  the  application  of  Theorem  2.  Bet  us  con¬ 
sider  the  following  matrix  characteristic  equation: 

+  B,s  +  C,  -  0  (11a) 


r*  .  mm _  _  DlttUHH 

D„  -  At  -  r  1,  D„  «  Bt  =  P5 71  13  71,  necessa 

|_0  lj  L 13 . 7  7.3J  conditi 
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III  Necessary  Conditions 

In  this  section  we  establish  some  necessary  conditions  for  the 
stability  of  multivariable  systems.  The  failure  to  satisfy  the 
necessary  conditions  for  stability  is  equivalent  to  the  sufficient 
conditions  for  the  instability  of  the  same  systems;  i.e., 
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Thaarara  3.  If  j  A/.|  j '  —  I,  2,  . . n  are  symmetric  such  that, 
there  exists  one  { .1/ ,  |  i  —  1,2,  n  which  is  negative  definite, 
negative  semi-definite,  or  indefinite,  then  the  system  in  equation 
(2)  is  unstable. 

Proef.  Suppose  the  system  is  asymptotically  stable  and  one 
of  j  A/,)  is  negative  definite,  negative  semi-definite,  or  indefinite. 
Since  the  stability  is  invariant  under  the  linear  transformation 
and  the  matrix  F  in  equation  (7)  is  a  stable  matrix.  Let  us  con¬ 
sider  the  following  equation: 

XF  +  FTX  -  -  Q  (12) 

where  Q  is  a  matrix  defined  in  equation  (9b).  By  Thoerem  4  in 
Bellman  191  Ip.  239)  and  the  theorems  in  Anderson  (10]  and 
Barnett  till  (p.  861,  we  know  that  equation  (12)  has  a  unique 
solution.  Since  Q  is  positive  semi-definite  and  rank  [Q|  =  rank 
[if]  —  m,  we  conclude  that  the  solution  X  of  equation  (12)  is 
also  positive  semidefinite.  Furthermore,  X  is  positive  definite 
if  the  pair  IF,  RT]  is  observable.  It  is  easy  to  verify  that  the 
matrix  P  which  was  defined  in  equation  (86)  satisfies  equation 
(12).  Therefore  X  —  P  is  positive  semidefinite  or  positive  def¬ 
inite.  This  implies  that  at  least  one  of  the  |  Af  ,|  is  positive  semi¬ 
definite  and  others  positive  definite  or  all  positive  definite.  This 
contradicts  our  assumption  that  one  of  the  ]  A/,)  is  negative 
definite,  negative  semidefinite,  or  indefinite.  Hence  the  system 
in  equation  (2)  is  unstable  if  one  of  the  { Af  <)  is  negative  definite, 
negative  semidefinite,  or  indefinite. 

To  show  an  application  of  Theorem  3,  consider  the  example 
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Al  ^  +  Bl  %  +  C'v  “  0  (l3a) 

aP  at 
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49J 
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Applying  equation  (5)  yields  the  matrix  Routh  array  and  M  s 


Hence,  the  system  is  unstable. 

The  next  criteria  is  another  necessary  oondition  which  we 
state  as  follows. 

Theorem  *.  If  det  Dw  >  0  and  det  D,.„  <  0,  or  det  Du  <0 
and  det  £>,.»  >  0,  and  Dn,  D,»  are  defined  in  equation  (14),  then 
the  system  in  equation  (2)  is  unstable. 

Proof.  Since  the  system  in  equation  (2)  has  the  matrix  char¬ 
acteristic  equation  (/>(«)]  in  equation  (14),  then  we  expand  the 
det  [/)(«)].  We  find  the  constant  term  is  equal  to  det  B,  «  det 
D.  „.  If  det  Du  >  0  and  det  »  <  0,  this  implies  that  the  coef¬ 
ficient  of  the  polynomial  det  {£>(«)}  has  a  negative  sign.  We  can 
then  conclude  that  the  det  [D(s)J  =  0  has  a  solution  with  a 
positive  real  part.  Hence  the  system  is  unstable. 

IV  Conclusion 

Some  necessary  and  some  sufficient  conditions  have  been  de¬ 
veloped  for  the  stability  of  a  class  of  multivariable  systems.  A 
linear  block  transformation  has  been  derived  for  transforming 
the  coordinates  of  an  observable  block  companion  form  to  the 
coordinates  of  a  block  Schwarz  form.  The  new  method  has 
partially  extended  the  scalar  Routh  criterion  to  the  matrix  Routh 
criterion  to  a  class  of  multivariable  systems. 
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